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IV. Interaction with Other Federal Agencies: 

1. AFOSR: A report has been submitted to Major Daniel Fant, manager of the unsteady 
program at the AFOSR on September 7, 1993. The report summarizes the important 
research results in supersonic vortex breakdown on delta wings and vertical-tail buffet 
problem. It also covers the concepts and ideas that benefit air force applications. It 
covers the period of October 1 -September 30, 1993. Major Fant is supporting part of the 
present work for the FY 1994 at $35 K. 

2. Numerical Aerodynamics Simulation: 

A proposal for computing resources has been submitted to the NAS facilities at NASA 
Ames in November 1992. The proposal has been approved for 1,100 hrs of computational 
time. Additional 500 hrs have been added to our account in July 1993. A summary has 
been submitted in July 1993 for publication in the NAS Annual Summaries. 

3. The Principal Investigator helped researchers at two branches at NASA Langley and a 
researcher at the Air Force Academy in Colorado to correct boundary conditions for 
unsteady flow calculation in the computer code “CFL3D”. 

V. Graduate Students: 

The Principal Investigator has three graduate students who are working on their M.S. thesis and 
Ph.D. dissertation. They are supported under this grant with the funds from NASA Langley 
Research Center and AFOSR. These students are: 

V.l. Mr. Mark W. Flanagan (US Citizen): He is supported under the present grant and a 
supplement from the Aerospace department. He is expected to defend his M.S. thesis 
in October 1993. He thesis is titled “Computational Simulation of Vertical Tail Buffet 
in a Configured Duct.” He is staying for his Ph.D. degree in the Aerospace Engineering 
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Department. His Ph.D. dissertation will focus on simulation and optimization control 
of three-dimensional tail buffet phenomenon. 

V.2. Mr. Steven J. Massey (US Citizen): He is supported under the present grant and a 
supplement from the Aerospace department. He is expected to defend his M.S. thesis in 
March, 1994. His thesis is titled “Computational Simulation of Tail Buffet Using a Delta 
Wing- Vertical Tail Configuration.” He is staying for his Ph.D. dissertation which will 
focus on simulation and optimizational control of tail buffet for bending and torsional 
responses. 

V.3. Ms. Margaret Kalisch (US Citizen): She joined our research group in May 1993. 
She has a M.S. degree from the Naval Postgraduate School in Monterey. CA. She is 
working on her Ph.D. degree in Aerospace Engineering. She is currently supported 
under the present grant. Her Ph.D. dissertation will focus on multi-mode coupling 
(roll/pitch and/or yaw) of a delta wing during vortex breakdown in transonic regime. 
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Computational simulation of vertical tail buffet has been accomplished using a simple generic model. The 
model consists of a vortex-breakdown flow in a configured circular duct and a vertical flat plate which is placed 
as a cantilever downstream of the vortex-breakdown flow. Vortex-breakdown flow is generated through a shock 
wave at the duct inlet. The unsteady, compressible, full Navier-Stokes (NS) equations and the aeroelastic 
equation for bending vibrations are solved sequentially to obtain the unsteady vortex-breakdown loads on the 
plate and its aeroelastic deflections, respectively. The grid displacements are obtained using an interpolating 

equation. 

t. INTRODUCTION 


The ability of modem fighter aircraft to fly and 
maneuver at high angles of attack is of prime im- 
portance. This capability is achieved, for example 
in the F/A-18 Fighter, through the combination of 
the leading-edge extension (LEX) for a delta wing 
and the use of a single or twin vertical tail(s). How- 
ever, at some flight conditions, the vortices emanating 
from the highly-swept LEX breakdown before reach- 
ing the vertical tails which get bathed in a wake of 
highly unsteady, turbulent swirling flow. The vortex- 
breakdown flow produces unsteady loads which in 
turn produce severe buffeting of the vertical tails and 
have led to their premature fatigue failure. Exper- 
imental research work has been conducted on dif- 
ferent models of the F/A-18 to understand the flow 
physics, measure the tails’ response [1-4] and alle- 
viate the buffet phenomenon [5]. Cole, Moss and 
Doggett [6] have found that the buffet response of a 
1/6 size model of the F-18 airplane occurs in the first 
bending mode. 

The crucial point in predicting the buffet charac- 
teristics is the knowledge of the driving unsteady air 
loads associated with the vortex breakdown flow. Ed- 
wards [7] presented a good assessment of the compu- 
tational cost of this problem for a full fighter aircraft. 
The principal author of this paper has proposed two 
simple models to simulate and study the vertical-tail 
buffet problem at a substantially reduced computa- 
tional cost in comparison with the cost of solving for 
the flow around a full fighter aircraft. The basic con- 
cept behind these models is to be able to generate 
an unsteady, vortex-breakdown flow and to place a 


vertical tail, which is cantilevered, downstream of the 
vortex-breakdown flow. In this way, the buffet prob- 
lem is isolated from the whole configuration and the 
computational resources are focused on a small region 
for high resolution. The first proposed model con- 
sists of a configured duct in which the inlet swirling 
flow is forced to breakdown either through a shock 
wave [8] (for transonic and supersonic inlet flows) 
or through a gradual adverse pressure gradient that is 
generated by the duct wall (for subsonic inlet swirling 
flows). Downstream of the breakdown flow, a verti- 
cal cantilevered tail is placed. In the second model, 
the configuration consists of a delta wing at a criti- 
cal angle of attack which produces breakdown of the 
leading-edge vortex cores. Downstream of the break- 
down flow, a vertical single or twin-tail configuration, 
which is fixed as a cantilever, is placed. 

Computational simulation using the second 
model has been presented in Ref. [9] by Kandil, 
Kandil and Massey. In the present paper, the first 
model is used for the computational simulation of ver- 
tical tail buffet. The model consists of a configured 
circular duct with supersonic swirling flow which is 
forced to breakdown through a shock wave at the 
inlet. Downstream of the vortex-breakdown flow a 
plate is placed as a cantilever which is fixed at the 
duct wall. The initial flow conditions are obtained by 
solving the unsteady, compressible, full NS equations 
accurately in time while the plate is assumed to be 
rigid. Next, the NS equations are sequentially solved 
with the aeroelastic equation of bending vibrations 
for the plate deflection. The grid is displaced using 
an interpolation equation. 
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2. FORMULATION AND 
COMPUTATIONAL SCHEME 

The present multidisciplinary problem is solved 
using three sets of equations which consist of the 
unsteady, compressible, full NS equations; the aeroe- 
lastic equation for bending vibrations of a cantilever 
beam, which approximates the plate bending vibra- 
tions, and an interpolation equation to move the grid 
due to the plate deflection. The NS equations are 
given in Ref. [9] and are solved using an implicit, 
upwind, flux-difference splitting (Roe scheme), finite- 
volume scheme. The aeroeiastic equation for the 
bending vibrations is given by 


are used. More details of the solution are given in 
Ref. [9]. 


For the grid displacements, the following inter- 
polation equation is used for the x-coordinate which 
is measured from the plate surface 


-n + l _ T n 
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where X r or / is the maximum x coordinate (mea- 
sured from the plate surface) of a grid point at the 
right (r) or left (/) boundary of the computational 
domain. With Eq. (5), a point on the plate is dis- 
placed by the total plate deflection and a point at 
either boundary is not displaced. 
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where r is the radial distance from the fixed sup- 
port along the beam length / tl E the modulus of 
elasticity, 1 the area moment of inertia, m the mass 
per unit length and A r is the net surface pressure 
force per unit length. The dimensionless form of 
u\ r, E , I and m are given by 
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where r m d is the duct inlet radius, the freestream 
air density, a ^ the freestream speed of sound, d * the 
plate thickness, b* the plate width and the denotes 
dimensional quantities. The geometrical and natural 
boundary conditions are given by 


*o .., -£m = £<*.«>- &<M)-o 

(^) 

The solution to Eq. (1) is given by 


n 

u'(r,f) = (4) 

;=i 

substituting Eq. (4) into Eq. (1) and using the 
Galerkin method, the equation for the generalized co- 
ordinates, qj{t), is obtained. The solution for q : (t) 
is obtained either by using a convolution integral or 
by using a four-stage Runge-Kutta scheme. For the 
present problem, six natural modes of vibration, <pj. 


3. COMPUTATIONAL RESULTS 

A configured, circular duct which is similar to 
the one used by the senior author in Ref. [8] is em- 
ployed for the present problem. The duct length here 
is 3.9 dimensionless units. It is configured such that 
a supersonic inlet flow produces a shock at the in- 
let section which is strong enough to break down the 
inlet swirling flow. A half of the duct axial plane 
is shown in Fig. 1. As a cantilevered plate at the 
duct wall of length l t = 0.4 is placed at the axial 
station x = 2.0, multi-block grids are used. There 
are three blocks of grids; the first is to the left of 
the plate (0 < x < 2.0, r > 0.6) which consists 
of 182x48 grid points in the axial and radial direc- 
tions, respectively; the second is to the right of the 
plate (2.0 < x < 3.9, r > 0.6) which consists of 
120x48 grid points in the axial and radial directions, 
respectively; the third extends from the center line to 
the level of plate tip (0 < r < 0.6, 0 < x < 3.9) 
which consists of 302x64 grid points in the axial and 
radial directions, respectively. The minimum radial 
grid size at the center line is Ar mm = 3.5 x 10~ 3 
and the minimum axial grid size at the plate surface 
is Ax min = 3 x 10" 4 . Figure 1 shows the grid used. 

The problem is solved for a quasi-axisymmetric 
flow by forcing the components of the flowfield vec- 
tor to be equal on two meridian planes in close prox- 
imity of each other. The inlet flow Mach number is 
1.75 and the Reynolds number is 10 5 (based on the 
duct inlet radius). The inlet profile for the tangential 
velocity is given by 
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where C/oo = 174, r m = 0.2 and k< = 0.1. The 
maximum swirling ratio 0, is at r = 0.224 
and its value is kept at 0.32. The radial velocity, 
v , at the inlet is set equal to zero and the radial 
momentum equation is integrated to obtain the inlet 
pressure profile. Finally, the density p is obtained 
from the definition of the speed of sound for the 
inlet flow. With this compatible set of profiles, the 
computations are carried out accurately in time with 
At = 0.001 and without inserting the plate until t = 9 
dimensionless time units. 

The boundary conditions for this step consists 
of the inflow Riemann- invariant conditions, solid 
boundary conditions at the duct wall, symmetry con- 
ditions at the centerline and extrapolating the flow 
variables at the duct exit boundary. Figure 2 shows 
two snapshots of the streamlines and Mach contours 
at t = 3 and t = 8. It is noticed that an oscillating 
shock developed at the duct inlet (see Mach contours) 
and the unsteady vortex-breakdown bubbles devel- 
oped behind the shock. The present solution is the 
same as that obtained earlier with one block of grid. 

3.1 initial Conditions 

For / > 9, a flat plate is impulsively inserted at 
x - 2.0 and is cantilevered at the duct wall. The 
computations are carried out accurately in time with 
A/ = 0.001 and with the plate being considered as 
a rigid boundary. Solid boundary conditions are en- 
forced on the plate surfaces. Figure 3 shows snap- 
shots of the streamlines and the Mach contours at 
t - 10, 11, 17, 23, 29 and 32. The streamlines 
show the continuous evolution, convection, merg- 
ing and shedding of vortex-breakdown bubbles be- 
hind the oscillating inlet shock. They also show an 
unsteady, separated vortical region behind the plate. 
The Mach-contours figures show the oscillating shock 
at the inlet and unsteady shocks between the plate tip 
and the duct wall. Having produced a strong vortex- 
breakdown flow ahead and behind the plate, this so- 
lution serves as initial conditions for the next step. 
The next step is to couple the computational fluid 
dynamics with the aeroelastic equation for bending 
vibrations of the plate along with Eq. (5) for the grid 
displacements. 

3.2 Vortex-Breakdown Flow and Plate Deflection 

For t > 34, the plate is allowed to deflect in 


bending by solving Eq. (1) with the force per unit 
length, N{r,t\ obtained from the net surface pres- 
sure on the plate. Since the plate is oscillating, the 
kinematical boundary conditions are based on the rel- 
ative velocity of the fluid with respect to the plate. 
The dynamical boundary condition on the plate 
is no longer equal to zero. The pressure gradient is 
modified to j£ = -pd p • n, where a p is the plate 
acceleration which is equal to \ plate and ” * s the 
unit normal to the plate. The computations are carried 
out accurately in time with A t = 0.0001. The plate 
dimensions are l t — 0.4, d = 0.02 and 6 = 0.157. 
Its modulus of elasticity and its mass per unit length 
are E = 0.499 x 10 3 and m = 7.177, respectively. 
Figure 4 shows snapshots of the streamlines and 
Mach contours during one cycle of periodic response 
of the plate. The snapshots are marked by the num- 
bers 1, 2, 3 and 4 which correspond to the instants 1, 
2, 3 and 4 shown on Fig. 5. Figure 5 shows the force 
distribution on the plate, N % at different time levels 
(from t = 35 to t = 35.5). It also shows the plate 
deflections at its midpoint and its tip. It is noticed 
that the period of oscillation is 0.3 time units which 
corresponds to a frequency of 20.94. The maximum 
tip deflection of the plate is 2.5 x lO^ 4 , which is of the 
same order as that of the axial minimum spacing of 
the grid at the plate surface. 

4. CONCLUDING REMARKS 

A simple generic model is introduced to sim- 
ulate the vertical tail buffet problem due to vortex- 
breakdown flows. The model consists of a config- 
ured circular duct which produces unsteady vortex- 
breakdown flow for the inlet supersonic swirling flow 
due to an oscillating shock. Downstream of the 
vortex-breakdown flow, a plate is placed as a can- 
tilever at the duct wall. The problem is solved us- 
ing three sets of equations which include the NS 
equations, the aeroelastic equation for bending vi- 
brations and an equation for the grid displacements. 
The equations are solved sequentially using time ac- 
curate stepping to obtain the force distribution on 
the plate, the plate deflections and the grid displace- 
ments, respectively. It has been shown that the plate 
reaches periodic response in bending in a short time. 
The present work represents a proof of concept for 
a generic model to simulate vertical tail buffet. The 
next step is to consider the three-dimensional flow 
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problem with the plate surface placed along the ax- 
ial flow direction where both bending and torsional 
modes can be considered. 
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ABSTRACT 

Computational simulation of the vertical tail buf- 
fet problem is accomplished using a delta wing- 
vertical tail configuration. Flow conditions are se- 
lected such that the wing primary-vortex cores expe- 
rience vortex breakdown and the resulting flow in- 
teracts with the vertical tail. This multidisciplinary 
problem is solved successively using three sets of 
equations for the fluid flow, aeroelastic deflections 
and grid displacements. For the fluid dynamics part, 
the unsteady, compressible, full Navier-Stokes equa- 
tions are solved accurately in time using an im- 
plicit, upwind, flux-difference splitting, finite-volume 
scheme. For the aeroelastic part, the aeroelastic equa- 
tion for bending vibrations is solved accurately in 
time using the Galerkin method and a convolution 
integral solution. The grid for the fluid dynamics 
computations is updated every few time steps using 
an interpolation equation. The computational appli- 
cations include a delta wing of aspect ratio 1 and two 
cases for the vertical tail. The first case is for a tail 
of aspect ratio of 2 located 0.5 chord length from the 
wing trailing edge. The second case is for a tail of 
aspect ratio of 1 located at the wing trailing edge. 

INTRODUCTION 

Recently, the design of modem fighter aircraft has 
been focused on high angle of attack maneuverabil- 
ity at high loading conditions, renewing the interest 
in the tail buffet problem. For these fighters, the abil- 
ity to fly and maneuver at high angles of attack is of 
prime importance. This capability is achieved, for ex- 
ample in the F/A-l 8 fighter, through the combination 
of a leading-edge extension (LEX) and a delta wing 
and the use of vertical tails. The LEX maintains lift at 
high angles of attack by generating a pair of vortices 
that trail aft over the top of the aircraft. The vortex 
entrains air over the vertical tails to maintain stability 
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of the aircraft. This combination of LEX and verti- 
cal tails leads to the aircraft excellent high angle of 
attack performance. However, at some flight condi- 
tions, the vortices emananting from the highly-swept 
LEX of the delta wing breakdown before reaching the 
vertical tails which get bathed in a wake of highly- 
turbulent, swirling flow. The vortex-breakdown flow 
produces severe buffet on the vertical tails and has 
led to their premature fatigue failure. Buffeting of 
the vertical tails of the F/A-l 8 fits into this category. 
During flight operations of this airplane large vibra- 
tions of the vertical tails have been observed. 

Sellers et al. 1 conducted some three-component 
velocity surveys for a YF-17 model (a configura- 
tion similar to the F-18) at low speeds. Their re- 
sults clearly show that at 25° angle of attack the 
vortex produced by the LEX experiences breakdown 
and that there are large fluctuations in the velocity in 
the vicinity of the vertical tails. They measured rms 
fluctuations as high as 40% of the freestream stream 
velocity. Erickson, et al. 2 presented a wind tunnel 
investigation of the F/A-l 8 aircraft. The investiga- 
tion focused mainly on the measurements of steady 
forces and pressures on the LEX and laser light sheet 
measurements of the vortex structure. Some water 
tunnel studies conducted by Wentz^ using an F-18 
model also showed that the vortex produced by the 
LEX of the wing breaks down ahead of the vertical 
tails at angles of attack of 25° and higher. If these 
flows contain substantial energy at frequencies corre- 
sponding to the lower modes of vibration of the tail 
structure, significant structural response can result. 

Another wind tunnel investigation of buffeting is 
published by Lee and Brown 4 . The buffeting of the 
vertical fin of a rigid 6% model of the F/A-l 8 has 
been investigated. Unsteady pressure measurements 
on the vertical fin were conducted and the vortex flow 
structure behind the fin was studied. The study was 
carried out with LEX fences on and off to conclude 
that the LEX fence has a small influence on the 
steady balance measurements such as lift and pitching 
moment. 
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Rao, Puram and Shah 5 proposed two aerodynamic 
concepts for alleviating high-alpha tail buffet charac- 
teristics of LEX vortex dominated twin tail fighter 
configurations. These concepts were explored in low 
speed tunnel tests on generic models via flow visu- 
alizations, 6-component balance measurements and 
monitoring of tail dynamics. Passive dorsal-fin exten- 
sions of the vertical tails, and an active LEX arrange- 
ment with up-deflected edge sections were evaluated 
as independent means of re-structuring the adverse 
vortical flow environment in the tail region. Each 
of these techniques successfully reduced the buffet. 
Used in combination, the two concepts indicated sig- 
nificant tail buffet relief with relatively minor impact 
on the high a configuration aerodynamics. 

Cole, Moss and Doggett 6 tested a rigid, 1/6 size, 
full-span model of an F-18 airplane that was fit- 
ted with flexible vertical tails of two different stiff- 
ness. Vertical-tail buffet response results were ob- 
tained over the range of angles of attack from —10° 
to +40°, and over the range of Mach numbers from 
0.3 to 0.95. Their results indicated that the buffet 
response occurs in the first bending mode, increases 
with increasing dynamic pressure and is larger at M 
= 0.3 than that at a higher Mach number. 

So far, there are no available theoretical or com- 
putational methods to predict and control the aeroe- 
lastic buffet characteristics of vertical tails. The cru- 
cial point in predicting the buffet characteristics is 
the knowledge of the driving unsteady airloads asso- 
ciated with flow separations and vortex breakdown. 
Edwards 7 presented a good assessment of the compu- 
tational cost of this problem for a full fighter aircraft. 
The principal author of this paper has proposed two 
simple models to simulate and study the vertical-tail 
buffet problem at a substantially reduced computa- 
tional cost in comparison with the cost of solving for 
the flow around a full fighter aircraft. The basic con- 
cept behind these models is to be able to generate 
an unsteady, vortex-breakdown flow and to place a 
cantilevered vertical tail, downstream of the vortex- 
breakdown flow. In this way, the buffet problem is 
isolated from the whole configuration and the com- 
putational resources are-focused-on a small region for 
high resolution. The first proposed model consists of 
a configured duct in which the inlet swirling flow is 
forced to breakdown either through a shock wave 8 
(for transonic and supersonic inlet flows) or through 
a gradual adverse pressure gradient that is generated 
by the duct wall (for subsonic inlet swirling flows). 


Downstream of the breakdown flow, a vertical can- 
tilevered tail is placed. In the second model, the con- 
figuration consists of a delta wing at a critical angle 
of attack which produces breakdown of the leading- 
edge vortex cores 9 . Downstream of the breakdown 
flow, a vertical single or twin-tail configuration which 
is fixed as a cantilever is placed. 

In the present paper, the second model is con- 
sidered for the computational simulation. The model 
consists of a delta wing of aspect ratio 1 and two cases 
of a single vertical tail. In the first case, the tail is of 
aspect ratio 2 and is placed at 0.5 root-chord length 
from the wing trailing edge. In the second case, the 
tail is of aspect ratio 1 and is placed at the wing trail- 
ing edge. The tail is cantilevered at the lower side and 
the configuration is pitched at a 35° angle of attack. 
The flow Mach number and Reynolds number are 0.4 
and 10,000, respectively. The problem is solved se- 
quentially using the time-accurate integration of the 
laminar, unsteady, compressible Navier-Stokes equa- 
tions, the aeroelastic equation for a beam in a bending 
mode and an equation to update the locations of grid 
points. 

FORMULATION 

The formulation of the problem consists of three 
sets of governing equations along with certain ini- 
tial and boundary conditions. The first set is the 
laminar, unsteady, compressible, full Navier-Stokes 
equations. The second set consists of the aeroelas- 
tic equations for the vibration modes which could be 
coupled bending and torsion modes. In the present 
paper, only the bending vibration is considered with- 
out structural damping or nonlinearities, which could 
be added in the future without any difficulty. The 
third set consists of an equation for calculating the 
grid displacements due to the tail deflections. The 
literature shows various methods for computing the 
grid displacements. The simplest methods use sim- 
ple interpolation functions such that the grid points 
adjacent to the aeroelastic surface move with the sur- 
face while the grid points at the computational-region 
boundary do not move 10 ' 11 . In the more advanced 
methods for moving the grid, the grid is simulated 
as a static 12-14 or dynamic truss. The unsteady, 
linearized, Navier-displacement equations have also 
been used successfully by Kandil et al. to move the 
grid dynamically 15 ' 16 . In the present paper, we use 
simple grid interpolation to move the grid. Next, the 
governing equations for each set are given: 
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Fluid Flow: 

The conservative form of the dimensionless, 
unsteady, compressible, full NS equations in 
terms of time-dependent, body-conformed coordi- 
nates f 1 ,^ 2 and f 3 is given by 

dQ'dErr, 3, 5 =1-3(1) 


— -+■ 
dt d£ 

where 


The reference parameters for the dimensionless form 
of the equations are c*, floe - f * / Q-x « Poe and p, x for 
the length, velocity, time, density and molecular vis- 
cosity, respectively. The Reynolds number is defined 
as Re = PocVoccVMx., where c* is root-chord length 
of the wing. The pressure, p, is related to the total 
energy per unit mass and density by the gas equation 


d£ s 


P = ( 7 - 1 )P 


Q = j = j[p,pui,pw2,pu 3 ,pe]' 

E m = inviscid flux 

= j[pUm,PU\U m + d\Z m p,pU2Um 
+ d2t m p,pu 3 U m + d 3 t m p,{pe + p)Un 

dCa 

-~dT p] 


( 2 ) 

(3) 


€ ~ 2 ( u i + u 2 + u i) 


(9) 


(4) 


(E v ) s = viscous and heat-conduction flux in £ s 
direction 

= j[0 ,d k Z s T kl ,dkt s T k 2,d k Z s TM, 

d k e(u n T kn -q k )} 1 ; k= 1-3, n= 1-3 (5) 

dC 


Um = d k C*k + 


dt 


( 6 ) 


The first element of the three momentum elements of 
Eq. (5) is given by 

s _ AW 
8k * Tkl = ~rT 


The second and third elements of the momentum 
elements are obtained by replacing the subscript 1, 
everywhere in Eq. (7), with 2 and 3, respectively. 
The last element of Eq. (5) is given by 

MooP 


d k t, s { u p T kp Qk) — j^ e 

- l d pt Sd kZ n ) u p-{j^ 
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(7 “ l ) Pr 




,d(a 2 ) 1 


dC 


The viscosity is calculated from the Sutherland law 

" = r3/2 (H£) ,c = 0,4317 (l0) 

and the Prandtl number P r = 0.72. In Eqs. OM8), 
the indicial notation is used for convenience. 

Aeroelastic Deflections: 

In the present paper, the vertical rectangular tail 
is treated as a homogeneous, uniform, cantilevered 
beam with rectangular section. For bending vibration, 
the dimensionless equation for the deflection w(:,t ) 
is given by 


d 4 w d 2 w . 

EI a? + m W = N{z ' t] 


(H) 



w = 

d k ediC - Idtf'dkC 

\ du k 

I dC I = 

dull 

(7) where c* 


where z is the vertical distance from the fixed sup- 
port along the beam length, /<, E the modulus of 
elasticity, / the area moment of inertia, m the mass 
per unit length and N is the net surface pressure 
force per unit length. The dimensionless form of 
w, 2 , E , /, m and N are given by 

w* _ z* r _ £* 

~c* , 2 - c*’ p* c a£' 

1 = 1 a. = Jr (i2) 

Poo a oc c 
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sound, d * the tail thickness, b* the tail chord length, 
A r * the net surface pressure force per unit length and 
the denotes dimensional quantities. The geomet- 
rical and natural boundary conditions are given by 

dw , , d 2 w . ,, _ 0V. _ n 

t«(o,t) - -r— (o,<) — -tty ( h,t) — 0,3 (^’0 ~ ® 
dz dz d, (13) 


;p = 1 - 3 


( 8 ) 


The solution to Eq. (11) is given by 

n 

w(z,t) = ( i4 ) 

j=i 
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where <pj(z ) are comparison functions which satisfy 
the boundary conditions, and they are given by the 
natural modes of vibration 


4>j{z) 


sin 0jZ - sinh/? r 2 + 
sin 0jlt + sinh/? r ff 


cosh QjZ — cos 0jZ 
cosh 3 , 1 1 + cos 3j!j 

(15) 


and 0j are the eigenvalues obtained from the solution 
of cos 0li cosh 0li = -1. Substituting Eq. (14) into 
Eq. (11) and using the Galerkin method, the follow- 
ing equation is obtained for the generalized coordi- 
nates qj(t): 


n n 

y m rjqj 3- y ^ k r jqj = N r (t) ; r 
j - i i =1 

where 


1,2 ,... ,n 
(16) 


m r j = elements of mass matrix 


= m I 4> r 4>j 

Jo 


dz 


(17.a) 


k r j = elements of sti f fness matrix 


= EI j 

Jo 


d 2 <p r d 2 <t>j 
dz 1 dz 2 


dz 


(17. b) 


t 1 ' 

N r = generalized force = I N<f> r dz (17. c) 

Jo 

Since m, I and E are constants and <t> r are orthogo- 
nal, the mass and stiffness matrices are diagonal and 
the set given by Eq. (16) is decoupled. Hence, the 
solution to the present simple case reduces to the so- 
lution of a decoupled set of second-order, ordinary- 
differential equations, where each equation corre- 
sponds to one of the natural modes. The solution 
is obtained either by using a closed-form convolution 
integral or by using a four-stage Runge-Kutta scheme. 
For the jth mode shape, Eq. (16) becomes 


, ^rr 

Qj 4 Qj — 

m rr 


Nr(t) 
m rr 


(no summation over r) 


(17) 


The convolution integral solution for this equation is 
given by 


1 f 1 - 

qj(t) = / A' r (C)sinuv(* - Od( 

uj r m rr J 0 

q,(0) 

+ qj(0)cosLj r t H sin uj r t (18) 

U>r 


where u: 2 = ^7 = ^rS9j(0) is the initial displace- 
ment and qj( 0) is the initial velocity. 


If m, / and E are functions of z, then Eq. (16) 
will be a coupled set and one needs to use a normal- 
mode shape transformation from q(t) to rj(t) to de- 
couple the resulting set. For a coupled bending and 
torsion vibration, the present procedure can be gen- 
eralized to obtain the solution. For the aeroelastic 
equations in the latter case, the solution is obtained 
using the four-stage Runge-Kutta scheme. 


Grid Displacements: 

In the present paper, we use simple interpolation 
functions to displace the grid due to the tail deflection. 
For bending vibrations, the tail deflection at any point 
on its surface, is in the spanwise direction (y 

coordinate). The spanwise coordinate of a grid point 
at the time level n + 1 is computed from the 
equation 


(,9) 


where Y i Jdim / 4 t is the maximum y coordinate of 
a grid point at the boundary of the computational 
domain with an index of JdimjA. Thus, the tail 
displacement is distributed through a cosine 

function among the y coordinates of the grid. Thus, 
a point on the tail is displaced by the total deflection 
and a point at the boundary is not displaced. 


Boundary and Initial Conditions: 

Boundary conditions consist of conditions for the 
fluid flow and conditions for the aeroelastic deflection 
of the tail. For the fluid flow, the Riemann-invariant 
boundary conditions are enforced at the inflow and 
outflow boundaries of the computational domain. At 
the plane of geometric symmetry, a periodic boundary 
condition is specified with the exception of grid points 
on the tail. On the wing surfaces, the no-slip and 
no-penetration conditions are enforced and = 0. 
On the tail surfaces, the no-slip and no-penetration 
conditions for the relative velocity components are 
enforced (points on the tail surface are moving). The 
normal pressure gradient is no longer equal to zero 
since points on the tail surface are accelerating. This 
condition becomes ^ = — Poo^t • Ti . where d< is the 
acceleration of a point on the tail, which is equal 
to For the boundary conditions of the 

aeroelastic deflection of the tail, they are given by 
Eq. (13). 

Initial conditions consist of conditions for the fluid 
flow and conditions for the aeroelastic deflection of 
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the tail. For the fluid flow, the initial conditions 
correspond to the freestream conditions with u\ = 
U 2 = U3 = 0 on the wing and tail. For the aeroelastic 
deflection of the tail, the initial conditions for any 
point on the tail is that the displacement and velocity 
are zeros, w(z,o) = ^(2,0) = 0 

METHOD OF SOLUTION 

The first step is to solve for the fluid-flow prob- 
lem using the vortex-breakdown conditions and keep- 
ing the tail as a rigid beam. Equations (1)-(I0) are 
solved time-accurately using the implicit, upwind, 
flux-difference splitting finite- volume scheme. The 
grid speed is set equal to zero in this step. This 
step provides the flow field solution along with the 
pressure difference across the tail. The pressure dif- 
ference is used to generate the force per unit length of 
the tail, N. Next, Eqs. (14) and (18) are used to ob- 
tain the tail deflections, Equation (19) is used 

to compute the grid coordinates. The metric coeffi- 
cients of the coordinate Jacobian matrix are updated 

Arm 

as well as the grid speed, Next, the compu- 

tational cycle is repeated using Eqs. (1)-(10) for the 
pressure difference across the tail, Eqs. (14) and (18) 
for the tail deflections and Eq. (19) for the grid co- 
ordinates. It should be noted that the time step for 
the fluid-flow problem, At f, is an order of magnitude 
less than the time step for the aeroelastic deflection, 
At,j. Moreover, the maximum tail deflection tv for 
each A is very small. Hence, one does not need 
to compute w for each time step, A tj. For example, 
if A t,j = 10 A tj, the computation of the aeroelas- 
tic deflections and grid coordinates can be performed 
every 10 A tj. 

COMPUTATIONAL APPLICATIONS 
CASE 1: 

The delta wing-vertical tail configuration consists 
of a sharp-edged, delta wing of aspect ratio 1 and 
a rectangular, vertical tail of aspect ratio 2, which 
is placed in the plane of geometric symmetry. The 
vertical-tail leading edge isJocated .at .0.5 roo(-chord 
length from the wing trailing edge. The lower edge 
of the tail is along the wing axis and the tail is 
clamped at that edge. The wing angle of attack is 
35° and the freestream Mach number and Reynolds 
number are 0.4 and 10,000, respectively. An O-H 
grid of 125x85x84 grid points in the wrap-around. 


normal and axial directions, respectively, is used for 
the solution of the fluid-flow part of the problem. 

The solution and analysis of this problem have 
progressed through rigorous steps in order to prove 
the capability of the present model for simulating the 
present multidisciplinary- problem: 

Step 1. Fluid Flow Around the Configuration 
(Initial Conditions): 

The laminar, unsteady, compressible, full Navier- 
Stokes equations have been integrated time accurately 
using the implicit, upwind, flux-difference splitting, 
finite-volume scheme with At = 0.003. During this 
step, the tail is kept rigid. The results of this step are 
used as initial conditions for the second step. Fig- 
ure 1 shows the spanwise, surface-pressure coefficient 
(C p ) on the wing at different chord stations at it = 
12,400 (37.2 dimensionless time). At x = 0.421, the 
Cp-curve shows asymmetry with the pressure on the 
left side having less suction than the pressure on the 
right side (looking in the upstream direction). This 
indicates that the left-vortex-core size is enlarging. In 
Fig. 4, it is noticed that at this location a spiral saddle 
point is formed and a spiral vortex-breakdown mode 
is developing downstream of that point. Returning 
back to Fig. 1, it is noticed that the asymmetry ex- 
ists in all the cross-flow planes downstream of x = 
0.421. At x = 0.814, a rapid decrease in the suction 
pressure on the left side is observed and Fig. 4 shows 
several spiral saddle points at that location. At x — 
0.944, a rapid decrease in the suction pressure on the 
right side is observed and spiral vortex breakdown 
develops in the right vortex core (see Fig. 4). Figure 
2 shows the total-pressure-loss contours and instan- 
taneous streamlines for the left and right sides at dif- 
ferent chord stations. The asymmetry is very clear at 
x = 0.65 and 1.0. At x = 1.0 and on the right side, it 
is noticed that a repelling spiral saddle point exists, 
which indicates the existence of vortex breakdown. 

Figure 3 shows the instantaneous streamlines in 
cross-flow planes at two chord stations x — 1 .375 and 
x = 2.019. The vertical-tail leading edge is located 
at x = 1.5 and its trailing edge is at x = 2.0. It is 
clear that the vertical tail is subjected to asymmetric, 
unsteady, surface-pressure loads due to the vortex- 
breakdown flow. Figures 4 and 5 show a plan view 
of the wing and a three-dimensional view of the 
wing-tail configuration, respectively, along with the 
spiral saddle points and the asymmetric spiral vortex- 
breakdown modes of the leading-edge vortex cores. 
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Figure 6 shows snapshots of the pressure contours 
on the right (J — 1) and left (J = Jdim) surfaces 
of the vertical tail and the pressure-difference on it 
at two time levels; it = 12,400 and it = 12,600. A 
close inspection of these loads reveals that the tail 
is subjected to both bending and torsional loads. In 
the present paper, only the bending vibration is taken 
into consideration (torsional vibrations are also under 
consideration). Figure 7 shows the corresponding 
lumped, aerodynamic force per unit length of the tail, 
N(z), at two time levels, / = 37.2 and t = 37.8. 
It is noticed that within this short time of change, 
N(z) changes rapidly in magnitude and direction. 
The results at it = 12,600 are used for the initial 
conditions of the buffet problem (i.e. letting the tail 
deflect and interact dynamically with the flow). 

Step 2. Fluid Flow/Tail Deflections Interaction: 

With the initial-flow conditions obtained in Step 1 
at it = 12,600, the aerodynamic force per unit length 
N r (z,t) is provided discretely to the aeroelastic pro- 
gram and the generalized forces N r are computed 
for six natural mode shapes. The deflections w are 
computed using Eqs. (14) and (18). Next, Eq. (19) 
is used to obtain the updated y-coordinates for the 
grid, the metric coefficients of the Jacobian matrix, 
the grid speed; and the velocity and acceleration of 
the points on the tail for the modified boundary condi- 
tions for the fluid-flow solution. Next, the fluid flow 
part is solved and N r (z,t ) is obtained and the com- 
putational cycle is repeated. The tail is treated as a 
homogeneous, uniform, rectangular beam with rect- 
angular cross section of thickness d = 0.01 and width 
b = 0.5. The dimensionless modulus of elasticity is 
E - 0.4903 x 10 2 and El = 0.0204 x 10 _4 . The 
dimensionless mass per unit length is m = 11 .428. 
Figure 8 shows the histories of tail deflection, w, and 
net pressure force per unit length, N , versus the ver- 
tical coordinate z every 500 time steps (every 1.5 
dimensionless time). The deflection and net pressure 
force versus time are also shown for the tail midpoint 
(z = 0.5) and its free-end point z — 1.0. The time- 
step counter n is for the interaction case only without 
the time steps for initial conditions (it = 12,600). The 
early time levels show very large N values with larger 
values at r = 0.5 than at z — 1.0. The tail deflection 
w shows a growing oscillatory response reaching a 
maximum value for the tail free end at n = 3,700 (it 
= 16,300). Then, the total deflection shows a damped 
oscillatory response due to the aerodynamic damping 


of the flow. The aerodynamic damping is automati- 
cally included in these results due to the method of 
solution which considers the fluid flow, the tail aeroe- 
lastic deflection and the grid deformation simultane- 
ously. The solutions do not show periodic responses 
until n = 6,500 (it = 19,100). It is also noticed that 
all the tail deflections are positive. 

Figure 9 shows cross-flow velocities and instan- 
taneous streamlines along the D-ray plane on the 
right and left sides showing the asymmetric vortex- 
breakdown flows of the leading-edge vortex cores at 
n = 6,500 (it = 19,100). In Fig. 10, the instanta- 
neous streamlines are shown in cross-flow planes be- 
fore and through the vertical tail at chord stations of 
x = 1.37, 1.51, 1.67, and 1.8. These sections show 
the asymmetric vortex-breakdown bubbles adjacent to 
the deflected tail. Figure 1 1 shows snapshots of the 
surface pressure contours at two instants, n = 3,150 
and 6,500 (it = 15,750 and 19,100) on the two sides 
of the vertical tail; right side is at J = 1 and left side 
is at J = Jdim. This figure shows the asymmetric 
unsteady pressure loads which affect the tail. 

To show that the tail deflection and its dynamic 
response affect the flowfield upstream of the tail loca- 
tion, the surface-pressure-coefficient distributions at 
several chord stations of the wing at it = 19,100 are 
shown in Fig. 12. Comparing the Cp’s of Fig. 1 (it 
= 12,400) and those of Fig. 12, it is observed that 
the Cp curves as of x = 0.421 are different. Since 
the flowfield is subsonic throughout, the disturbances 
created by the tail deflection and dynamics propa- 
gate upstream affecting the flowfield and the vortex- 
breakdown critical points. Figure 13 shows deformed 
grids in a cross-flow plane passing through the tail at 
x = 1.75 at it = 15,750 and 19,100. The grid defor- 
mations and the tail deflection are clearly noticed. 

In Fig. 14, color graphics of a top view, tilted 
front view and a three-dimensional view of the wing- 
tail configuration are shown at it = 19,100. The top 
view shows the locations of the spiral saddle points 
on the right and left sides of the wing (their locations 
are different from those of Fig. 4). The tilted front 
view shows the trace of the deflected tail and the 
total-pressure surfaces which show the breakdown of 
the leading-edge vortex cores. The three-dimensional 
view shows the whole configuration including the 
total pressure surfaces, the spiral saddle points and 
the surface-pressure contours (compare with those of 
Fig. 5). 
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CASE 2 

With the same delta wing and same flow condi- 
tions, the vertical tail aspect ratio is reduced to one 
(b = 0.5, h = 0.5) and its location is changed to the 
wing trailing edge from x = 1 .5 to x = 1 .0. The solu- 
tion for the initial conditions is repeated for this new 
configuration until it = 12,600. The fluid flow, tail 
deflection and grid displacement interactions are car- 
ried out for n = 4,500 (it = 17,100) with A / = 0.003. 
The material properties of the tail are kept fixed. 

Figure 1 5 shows the histories of tail deflection and 
net pressure force per unit length versus the vertical 
coordinate z every 500 time steps (every 1 .5 dimen- 
sionless time). The deflection and net pressure force 
versus time are also shown for the tail midpoint (~ = 
0.25) and its free-end point (z = 0.5). The tail deflec- 
tion shows a growing oscillatory response reaching 
a maximum value for the midpint at n = 1,200 (it 
= 13,800). Then, the tail deflection shows a rapidly 
damped oscillatory response due to the aerodynamic 
damping of the flow. The tail deflections lead the 
net pressure forces and the phase lead is more pro- 
nounced at the free end of the tail. The solutions do 
not show a periodic response until n = 4,500 (it = 
17,100). However, the tail deflections show positive 
and negative values. The tail deflections are also no- 
ticed to change from the second bending-mode shape 
to the first bending-mode shape. 

Figure 1 6 shows cross-flow velocities and instan- 
taneous streamlines along the D-ray plane on the 
right and left sides showing the asymmetric vortex- 
breakdown flows of the leading-edge vortex cores at 
n = 4,500 (it = 17,100). In Fig. 17, the instantaneous 
streamlines are shown in cross-flow planes passing 
through the tail at a - = 1 .0 1 , 1.18, 1 .38 and 1 .5. These 
sections show the asymmetric vortex-breakdown bub- 
bles adjacent to the tail surface and covering a larger 
height of the tail than those of Fig. 10 (the tail height 
here is one-half the tail height of Case 1). 

In Fig. 18, color graphics of a top view, a tilted 
front view and a three-dimensional view of the wing- 
tail configuration are shown at it = 17,100. Compar- 
ing these views with those-of Fig. 14,- it is observed 
that the spiral saddle points of vortex breakdown, sur- 
face pressure and total-pressure surfaces are different 
from those of Fig. 14. Again, this shows conclusively 
that location, shape, deflections, and dynamics of the 
tail substantially affect the flowfield upstream of the 
tail, i.e.; on the wing. In Fig. 8, the tail deflection 


is noticed in the top view and the three-dimensional 
view as well. 

CONCLUDING REMARKS 

The buffet problem of a vertical tail due to the 
interaction of vortex-breakdown flow with the tail 
has been simulated computationally and efficiently 
using a delta wing-vertical tail configuration. The 
wing aspect ratio and flow conditions have been care- 
fully selected in order to produce an unsteady vortex- 
breakdown flow. The solution has demonstrated the 
development of the tail buffet due to the unsteady 
loads produced by the vortex-breakdown flow. The 
problem is a multidisciplinary problem which re- 
quires three sets of equations to obtain its solution. 
The first set is the unsteady Navier-Stokes equation 
which is used to obtain the aerodynamics force per 
unit length on the tail. The second set is the aeroe- 
lastic equation for bending vibrations which is used 
to obtain the tail deflections, velocities and accelera- 
tions. The third set is the grid displacement equation 
which is used to update the grid coordinates due to 
the tail deflections. The three sets of equations are 
solved sequentially and accurately in time. The com- 
putational applications included two cases of delta 
wing-vertical tail configurations. Fixing the flow con- 
ditions and the geometry of the delta wing, two tail 
aspect ratios and locations are used. Initially, tail de- 
flections and aerodynamic loads were higher for the 
first case than the second case. Later on, the deflec- 
tions were damped due to the aerodynamic damping 
of the flow. The solutions show that the tail location, 
shape, deflections and dynamics affect the flowfield 
upstream of the tail. Work is underway to include tor- 
sional modes to the bending modes, upgrade the tail 
model from beam equations to plate equations, and 
consider the wing-twin vertical tail configuration. 
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Fig. 4. Top view showing surface-pressure contours and streamlines in the vortex 
cores showing their spiral-vortex breakdown, AR W = 1. ARf 2. 
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Fig. 5. Three-dimensional view of delta wing-vertical tail configuration showing 
the spiral vortex breakdown of vortex cores, ARw — 1, ARt = 2. 
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Fig. 7. Lumped force per unit length on the tail at two time levels; 
it = 12,400; it = 12.600. 
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Fig. 10. Cross-flow instantaneous streamlines at different chord stations before 
and through the tail, it = 19,100, At = 0.003, ARt = 2. 
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ABSTRACT 

Transonic flow over a 65-degree swept-back, sharp- 
edged, cropped delta wing is investigated computationally 
HQing the time-accurate solution of the unsteady, com- 
pressible, full Navier-Stokes equations with an implicit, 
upwind, flux-difference splitting, finite-volume scheme. 
Coarse and fine O-H grids are used to obtain the so- 
lution. The grid consists of 125 x 85 x 84 points in the 
wrap-around, normal and axial directions, respectively. 
The results are presented for an angle of attack of 20° , 
March number of 0.85 and Reynolds number of 3.23 x 10 6 
(based on the wing chord length). With the fine grid, the 
results show that a system of shocks has been captured 
over the upper wing surface, and that the leading-edge 
vortex core experiences an unsteady, supersonic vortex 
breakdown after passing through a spanwise shock (ter- 
minating shock) near the wing trailing edge. The com- 
puted results at a certain time are in good agreement with 
the experimental data. Topological aspects of the vortex- 
breakdown flowfield are also presented and discussed. 

INTRODUCTION 

At sufficiently high angles of attack, vortex break- 
down for incompressible flows around delta wings has 
been observed along the leading-edge primary vortex 
cores. Two distinct forms of vortex breakdown have been 
documented experimentally 1 . The first form is the bubble 
type and the second form is the spiral type. The bub- 
ble type shows an almost axisymmetric sudden swelling 
of the vortex core into a bubble, while the spiral type 
shows an asymmetric, spiral, vortex filament followed by 
a rapidly spreading turbulent flow. Both types are charac- 
terized by an axial stagnation point and a limited region 
of reversed axial flow. Much of our knowledge of in- 
compressible vortex breakdown has been obtained from 
experimental studies of pipe flows where both types of 
breakdown and other types as well were generated and 
documented 2 " 4 . 

The major effort of computational study of vor- 
tex breakdown flows has also been -focused on iso- 
lated swirling flows. For incompressible flows, quasi- 
axisymmetric, bubble-type, vortex-breakdown flows were 
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computed using the Navier-Stokes equa t ions^ . Three- 
dimensional bubble and spiral vortex-breakdown flows 
were also computed for isolated swirling flows using 
the three-dimensional Navier-Stokes equations in the 
vorticity-velocity form or the primitive variables form*" 11 . 

Interaction between a longitudinal vortex and a trans- 
verse shock wave occurs in several flow applications 
which include transonic and supersonic flows over a delta 
wing or a strake-wing configuration at moderate to high 
angles of attack, a supersonic inlet ingesting a vortex, 
and a supersonic combustor where fuel is injected in a 
swirling jet to enhance fuel-air mixing 13-14 . For delta 
wings and strake-wing configurations, vortex breakdown 
is an undesirable phenomenon since it produces wing 
stall. Therefore, its occurrence needs to be delayed with 
passive or active control methods in order to increase 
the wing performance at large angles of attack. For 
a supersonic combustor, vortex breakdown is desirable 
since it enhances mixing of air and fuel and stabilizes the 
flame 15,16 . Therefore its occurrence needs to be enhanced 
and controlled. 

For supersonic flows, quasi -axisymmetric bubble-type 
vortex-breakdown 17-19 and three-dimensional bubble-type 
and spiral-type vortex breakdown 20 for isolated swirling 
flows have been recently computed by the present au- 
thors. Using compatible, inlet boundary conditions, the 
time-ac cu rate solutions of the unsteady, compressible, full 
Navier-Stokes equations were obtained to study the ef- 
fects of Reynolds number, Mach number, swirl ratio, 
type of exit-boundary conditions and grid fineness and 
distribution on the vortex-breakdown modes for internal 
and external flows. Several modes of vortex breakdown 
which include transient single-bubble, transient multi- 
bubble, periodic multi-frequency multi-bubble, quasi- 
steady two-bubble cell and spiral-type vortex break- 
downs have been obtained 21 . For three-dimensional 
vortex-breakdown flows in a swirling, supersonic jet 
flow, topological aspects of the critical points in the 
vortex-breakdown region have been studied and com- 
pared with the available experimental incompressible 
vortex-breakdown topology. 

Recent experimental measurements 22 26 of transonic 



flows around a 65° swept-back, cropped delta wing show 
that shock wave formation is likely to occur underneath 
the leading-edge primary vortex core. In cross-flow 
planes perpendicular to the wing, the cross-flow beneath 
the primary vortex reaches supersonic speeds and a cross- 
flow shock develops beneath the primary vortex similar 
to the supersonic flow in a convergent-divergent nozzle. 
These measurements also show that a transverse shock 
“terminating shock" which might cause primary-vortex- 
core breakdown could develop in an analogous manner 
to the shock that terminates the two-dimensional super- 
sonic pocket on an airfoil. A complete reconstruction 
of the three-dimensional flow field on the delta wing in 
this region was not possible experimentally 22- 26 . Com- 
putational simulations for transonic delta-wing flows have 
been developed by using the Euler equations 27 ' 2 * and 
the thin-layer, Navier-Stokes equations 29 . The Euler- 
equations solutions were not capable of fully resolving 
the flow in the terminating shock region and the thin- 
layer, Navier-Stokes solution did not address that region. 

In the present paper, we consider the transonic flow 
around a 65° sharp-edged, cropped delta wing at an angle 
of attack of 20°, a Mach number of 0.85 and a Reynolds 
number of 3.23 x 10 s . The purpose of the present numer- 
ical simulation and study is to construct the flow field 
over the wing with particular emphasis of the vortex 
core-terminating shock interaction region. The laminar, 
unsteady, compressible, full Navier-Stokes equations are 
solved accurately in time with an implicit, flux-difference 
splitting, finite-volume scheme. The computations are 
carried out with time-accurate stepping on two O-H grids; 
a coarse grid and a fine grid. Both grids consist of 
125 x 85 x 84 points in the wrap-around, normal and axial 
directions, respectively. The main difference between the 
coarse and fine grids is the distribution of the grid points 
normal to the wing surface within the thin viscous layer 
(to be discussed later on). 

HIGHLIGHTS OF FORMULATION 

AND COMPUTATIONAL SCHEME 

The conservative form of the dimensionless, unsteady, 
compressible, full Navier-Stokes equations is used for the 
formulation of the problem. The equations are written in 
terms of the time-independent, body-conformed coordi- 
nates and £ 3 (Ref. 18). 

The implicit, upwind, flux-difference splitting, finite- 
volume scheme is used to solve the unsteady, compress- 
ible, full Navier-Stokes equations. The scheme uses the 
flux-difference splitting scheme of Roe which is based 
on the solution of the approximate one-dimensional, Rie- 
mann problem. In the Roe scheme, the inviscid flux dif- 
ference at the interface of computational cells is split into 
two pans; left and right flux differences. The splitting is 
accomplished according to the signs of the eigenvalues of 
the Roe averaged-Jacobian matrix of the inviscid fluxes 
at the cell interface. The smooth flux limiter is used to 
eliminate oscillations at locations of large flow gradients. 


The viscous- and heat-flux terms are linearized in time 
and the cross-derivative terms are neglected in the im- 
plicit operator and retained in the explicit terms. The vis- 
cous terms are differenced using a second-order accurate 
central differencing. The resulting difference equation is 
approximately factored and is solved in three sweeps in 
the and ( 3 directions. The computational scheme 
is coded in die computer program “FTNS3D" which is a 
modified version of the CFL3D-code. 

COMPUTATIONAL RESULTS 

A 65° swept-back, sharp-edged, cropped delta wing 
with zero thickness is considered for die computational 
solutions. The cropping ratio (dp length/root-chord 
length) is 0.15. The wing angle of attack is 20 s , and the 
freestream Mach number and Reynolds number (based 
on the root-chord length) are 0.85 and 3.23 x10 s , re- 
spectively. The reason behind the present, selected flow 
conditions is because of the uncertainty of the existing 
experimental data 22-26 about the structure of the down- 
stream flow field of the leading-edge vortex core. The 
experimental data shows that a supersonic flow region 
appears on the upper wing surface near the plane of sym- 
metry. This flow region is terminated by a transverse 
shock (known as a terminating shock) in a similar way to 
the shock that terminates a supersonic pocket on a super- 
critical airfoil 23 . 

Grid: 

An O-H grid of 125 x 85 x 84 in the wrap-around, 
normal and axial directions, respectively is used for the 
computational simulation. The computational domain 
extends two-chord length forward and five-chord length 
backward from the wing trailing edge. The radius of the 
computational domain is four-chord length. Two grids 
have been constructed using the same number of grid 
points. The first is called the coarse grid and the second 
is called the fine grid. For the coarse grid, the grid points 
in the cross flow planes have been distributed using a 
Joukowski transformation which produces a minimum 
grid size, normal to the wing surface, that varies from 
SxKT* at the leading edge to 3xlO" 2 at the plane of 
symmetry. For the fine grid, the elliptical grid lines in 
the cross-flow planes have been constructed such that the 
minimum grid size normal to the wing surface, stays 
constant at 5X10 -4 from the leading edge to the plane 
of symmetry. Figures 1 and 2 show three-dimensional 
shape of the coarse and fine grids and a cross-flow plane 
along with its blow-ups. 

Time-accurate integration of the laminar, unsteady, 
compressible, full Navier-Stokes equations has been car- 
ried out with At = 0.001 for the coarse grid and At = 
0.0002 for the fine grid. The results showed that the 
leading-edge vortex core passes through a terminating 
shock which causes the vortex core to breakdown. More- 
over, it is shown that the flow becomes unsteady behind 
the terminating shock. 
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Validation of Surface Pressure: 

Figure 3 shows a comparison of the computed, span- 
wise surface-pressure coefficient (Cp) at different chord 
stations for the fine and coarse grids with the experimen- 
tal data of Erickson 30 (R« = 3.23 xlO 6 ) and Hartmann 24 
(R« = 2.38 xlO 6 and 4.57 xlO 6 ). The computed results 
are selected at t = 3.6. Obviously, the coarse-grid C p - 
curves do not show the suction-pressure peak correspond- 
ing to the secondary vortex and the correct location of the 
suction-pressure peak corresponding to the primary vor- 
tex. The coarse-grid C p -curves are similar to those of the 
Euler-equations solution. Therefore, they are discarded 
in this paper. The fine-grid C p -curves show the correct 
location of the suction-pressure peak corresponding to the 
primary vortex and the suction-pressure peak correspond- 
ing to the secondary vortex. The fine-grid Cp-curves at 
x = 0.3, 0.6 and 0.8 are in fair to good agreement with 
the experimental data. For x = 0.9, the fine-grid C,,- 
curve shows a substantial, rapid increase in the pressure 
coefficient (a decrease in the suction pressure). Figure 
4 shows the total-Mach contours and the streamlines in 
cross-flow planes at the chord stations of x = 0.3, 0.6, 0.8, 
0.9, 0.97 and 1.0. At x = 0.3, 0.6 and 0.8, the total-Mach 
contours show an oblique shock under the primary vortex 
and a small subsonic region to the right of the shock. The 
streamlines show the secondary separation to the right of 
the shock. This separation is due to the shock interac- 
tion with the surface boundary-layer flow and is also due 
to the adverse, spanwise pressure gradient created by the 
primary vortex. At x = 0.9, the shock under the primary 
vortex becomes weak as observed in the total-Mach con- 
tours and the primary-vortex size increases. At x = 0.97, 
the shock under the primary vortex disappears and the 
primary vortex diffuses and reduces to a repelling focus 
as shown by the streamlines. At x = 1.0, the repelling 
focus becomes a repelling line. The details of the flow 
structure shown at x = 0.9, 0.97 and 1.0 indicate that 
the primary vortex is going through a breakdown mode 
which is caused by a transverse shock (terminating shock) 
between x = 0.8 and x = 0.9. 

Terminating Shock: 

To show that a terminating, transverse shock exists 
and has been captured computationally, the static-pressure 
contours and total-Mach-contours on two planes are com- 
puted and displayed in Fig. 5. In this figure, the static- 
pressure contours are shown on the wing surface and 
the plane of symmetry, and the total-Mach contours are 
shown on the third plane (k = 3) above the wing (in the 
viscous layer) and on the plane of symmetry. The plane 
of symmetry contours clearly show the location, shape 
and strength of the terminating shock. Moreover, the 
Mach contours show that a substantial supersonic pocket 
(bounded by the sonic line and terminating shock) ex- 
tends from the wing vertex to the shock location of x = 
0.83, which is in good agreement with the experimental 
data 30 , where the shock is located at x = 0.84. The com- 


puted results show that the shock is a normal shock with 
a height of 0.4 which is equal to one-half the wing span. 

In the spanwise direction, the shock foot print (shown on 
the Mach contours at k = 3) extends beyond the primary- 
vortex location. A A-type shape of the shock-system foot 
print, which on one side of the wing, consists of the ter- 
minating shock and the shock under the primary vortex 
that runs along a ray plane from the wing vertex, is seen 
on the Mach contours at k = 3. 

Figure 6 shows the position of the ray lines from the 
wing vertex (which are marked by the letters A-H) and 
the static-pressure curves along these lines. The static- 
pressure curves show the spanwise locations of several 
points on the foot-print line of the terminating shock. The 
terminating shock is clearly seen to run in the spanwise 
direction from the plane of symmetry to the wing leading 
edge. It reaches its highest strength from the location of 
the primary vortex to the wing leading edge (from line 
E to line H). 

Vortex-Breakdown Structure: 

Having established the shock system that consists of 
the shock under the primary vortex and the terminating 
shock, the focus is directed on the structure of the flow be- 
hind the terminating shock. In Fig. 7, we show the total- 
Mach contours and streamlines on a ray plane at the 0.658 
spanwise location, which passes through the leading-edge 
vortex core. Blow-ups of the velocity vectors and stream- 
lines on this vertical plane are also shown in Fig. 7. The 
streamlines figures clearly show a two-bubble cell vortex 
breakdown. This is a typical three-dimensional vortex- 
breakdown mode which consists of an attracting saddle 
point (front) a repelling saddle point (rear), an attracting 
focus (top) and a repelling focus (bottom). Such a break- 
down mode is similar to the one which was captured for 
an isolated supersonic vortex in an unbounded domain 
in Refs. 20 and 21. The location of the attracting sad- 
dle point is at 0.97 along the ray line, which corresponds 
to 0.87 along the axial direction. The attracting focus 
point is characterized with spiralling-in streamlines and 
the repelling focus point is characterized with spiralling- 
out streamlines. The Mach contours show that the front 
surface of the voitex-breakdown bubbles is enclosed by a 
hemi-spherical shape-like shock surface. Figures 12 and 
13 show details of the flow structure on the wing plan 
view, on the plane of symmetry and on the ray plane at 
the 0.658 spanwise location (marked as J = 16 on Fig. 13). 
These figures and discussion give a complete construction 
of the flow structure including the shock system and its 
interaction with the leading-edge vortex core which pro- 
duces vortex-breakdown of the two-bubble-cell mode. 

Unsteadiness of the Vortex-Breakdown: 

The computations have been carried out with time- 
accurate stepping beyond t = 3.6. Figures 8—1 1 show the 
results at t = 5.52. These results show that the terminating 
shock moves in the upstream direction and so is the 
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two-bubble-cell vortex breakdown behind the terminating 
shock. Figure 8 shows that the repelling focus is at x = 

0.88 instead of x = 0.97 (Fig. 4). Figure 9 shows that 
the terminating shock in the plane of symmetry is at x = 
0.685 instead ofx = 0.83 (Fig. 5). The shock decreases in 
height and its thickness increases. Figure 10 shows that 
the size of the two-bubble cell vortex-breakdown region 
increases in comparison with the size at t = 3.6 (Fig. 7). 
Upstream of the terminating shock the flow stayed steady 
without any change. 

Beyond the time t = 5.52, the upstream shock motion 
stopped and the motion reversed its direction to the down- 
stream. The computations were not carried out beyond 
this instant due to its impeding cost. The unsteadiness of 
the terminating shock and the vortex-breakdown region 
behind it have also been observed experimentally by Ban- 
nik and Houtmann 23 . They also observed that the flow 
upstream of the terminating shock stayed steady without 
any change. These experimental observations undoubt- 
edly support and validate our computational results. 

CONCLUDING REMARKS 

The laminar, unsteady, compressible, full Navier- 
Stokes equations are integrated time accurately using the 
implicit, upwind, flux-difference splitting, finite-volume 
scheme to study and construct the flow field structure of 
a transonic flow around a 65° sharp-edged, cropped-delta 
wing. A A-shock system, which consists of a ray shock 
under the primary vortex core and a transverse terminat- 
ing shock, has been captured. Behind the terminating 
shock, the leading-edge vortex core breaks down into a 
two-bubble cell type. The terminating shock and the vor- 
tex breakdown region behind it is time dependent and 
appears to be oscillatory. The flow field ahead of the ter- 
minating shock stays steady without any change. This is 
consistent with the fact that the supersonic pocket along 
with the terminating shock do not allow disturbances to 
propagate upstream. The present results have been vali- 
dated using the available experimental data and they are 
in good agreement The present paper gives a complete 
construction of the flow field over the wing surface and 
in particular the structure of the flow at the terminating 
shock and behind it. 
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coefficient at different chord stations; Moo = 0.85, a = 20°, t = 3.6 














Spin-ill* 


1 .00 



PRESSURE CONTOURS ON THE TT1NG SURFACE 



MACH CONTOURS ON A CONSTANT K PLANE 


0.80 


0.00 



0.00 

PRESS. CONTOURS ON THE PLANE OF SYMUTRY 



0.00 1 * 00 
MACH CONTOURS ON THE PLANE OF SYMMTRY 


Sonic Une 


Fig. 5. Static-pressure and Mach contours on the wing and the 
plane of symmetry; Moo = 0.85, a = 20°, t — 3.6 




Fig. 6. Ray lines on the wing surface and the static-pressure variation 
along them; Moo = 6.85, a = 20°, t — 3.6 
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Fig. 7. Total-Mach contours, streamlines and velocity vectors on a ray plane passing through 
the vortex-breakdown two-bubble cell; Moo = 0.85, o = 20 , t — 3.6 


Fig. 8. Total-Mach contours and streamlines in cross-flow 
plane; Moo = 0.85, a = 20°, t = 5.52 
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Fig. 9. Static-pressure contours on the wing and the plane of 
symmetry; Moo = 0.85, a = 20°, t = 5.52 
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Fig. 10. Ray lines on the wing surface and the static-pressure 
variation along them; Moo — 0.85, a = 20°, t = 5. 
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Fig. 11. Total-Mach contours, streamlines and velocity vectors on a ray plane passing 
through the vortex-breakdown bubbles; M 0 0 = 0.85, a = 20°, t = 5.52 






Fig. 12 Surface-pressure and Mach contours and particle trance on wing 
and symmetry planes; M 0 0 = 0.85, a - 20°, t = 3.6 
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ABSTRACT 

The effects of freestream Mach number and angle of 
attack on the leading-edge vortex breakdown due to the 
terminating shock on a 65-degree, sharp-edged, cropped 
delta wing are investigated computationally. The compu- 
tational investigation uses the time-accurate solution of 
the laminar, unsteady, compressible, full Navier-Stokes 
equations with the implicit, upwind, flux -difference split- 
ting, finite-volume scheme. A fine O-H grid consisting of 
125 x 85 x 84 points in the wrap-around, normal and axial 
directions, respectively, is used for all the flow cases. 
Keeping the Reynolds number fixed at 3.23x10®, the 
Mach number is varied from 0.85 to 0.9 and the angle 
of attack is varied from 20° to 24°. The results show that 
at 20° angle of attack, the increase of the Mach number 
from 0.85 to 0.9 results in moving the location of the ter- 
minating shock downstream. The results also show that 
at 0.85 Mach number, the increase of the angle of at- 
tack from 20° to 24° results in moving the location of 
the terminating shock upstream. The results are in good 
agreement with the experimental data. 


INTRODUCTION 

The literature shows that vortical flows around delta 
wings in the low-speed regime have received a substantial 
volume of experimental 1-4 and computational 5-9 research 
work. In the high angle of attack range, vortical flows in 
the low-speed regime are characterized with three types 
of boundary-layer separation, namely; primary , secondary 
and tertiary separations. The primary separated flow rolls 
up into a strong primary vortex core which produces a 
strong suction-pressure peak on the wing surface. The 
spanwise adverse-pressure gradient of the primary vortex 
causes the spanwise, outboard-moving, boundary-layer 
flow to separate forming a secondary vortex with opposite 
sense of rotation to and smaller strength than that of the 
primary vortex. The spanwise adverse-pressure gradient 
of the secondary vortex causes the spanwise, inboard- 
moving, boundary-layer flow to separate forming a ter- 
tiary vortex with same sense of rotation as and substan- 
tially small strength than that of the primary vortex. The 
spanwise surface-pressure curves are characterized with 
three suction-pressure peaks which varies in strength cor- 

■ Research AiiocUte, Dept, of Aenxpace Engineering. Member AIAA. 

~ Pmfeuor. Eminent SchoUr end Cheinnm of Dept of Aert»p«ce Engineer- 

ing. Associate Fellow AIAA. » 

*** Senior Research Scientist, Computational Aerodynamic Branch, Associate 

Fellow AIAA. 


Copyright C 1993 by Omm. Kondil. Publithed by Hie America 
SLfof Aeroruutjci md Anroruuuci, Inc. wrth penrnuion. 


responding to the locations of the primary, secondary and 
tertiary vortices. When the angle of attack reaches a crit- 
ical value, the axial-pressure gradient and the high swirl 
ratio of the primary vortex produce a stagnation point 
along the path line of the primary-vortex core, and vor- 
tex breakdown of the primary core develops. Depending 
on the swirl ratio, axial pressure gradient and Reynolds 
number, the primary-core vortex-breakdown mode might 
be a bubble type, a spiral type or a bubble-spiral type. 

As the freestream Mach number increases, the vortical 
flow around the delta wing changes substantially due to 
the compressibility effects. In the supersonic flow regime, 
shock waves appear beneath or above the primary vor- 
tex, depending on the freestream normal Mach number 
and normal angle of attack. Experimental data • and 
the computational results 12-14 have shown these types o 
vortical-flow structures. The foot print of these shock 
waves runs along a ray line from the wing vertex. If 
the shock wave is beneath the primary vortex, it interacts 
with the spanwise, outboard-moving, boundary-layer flow 
and causes, in addition to the adverse pressure gradient 
produced by the primary vortex, secondary-flow separa- 
tion. If the shock wave is above the primary vortex, it 
flattens the primary vortex and the spanwise surface pres- 
sure curve. Comparison of the surface pressure distribu- 
tion over a delta wing in low-speed and supersonic-speed 
regimes, shows that the suction-pressure peak correspond- 
ing to the primary vortex is lower for the supersonic flow 
than that for the low-speed flow. 

In the transonic-flow regime, research work on vor- 
tical flows around delta wings was given adequate at- 
tention only recently. Understanding the steady and un- 
steady, transonic, vortical-flow structures around delta 
wings in the moderate-high angle of attack range is im- 
portant for increasing the performance quality of the new 
generation of supermaneuver aircraft (e.g. YF22). Recent 
experimental measurements of transonic flows around a 
65° cropped delta wing 15-21 show that a complex shock- 
wave system appears over the upper wing surface. The 
shock-wave system consists of a ray shock wave beneath 
the leading-edge primary vortex and a transverse, time- 
dependent 16 , normal-shock wave (known as a terminat- 
ing shock) which runs from the plane of symmetry to the 
wing leading edge. The terminating shock wave interacts 
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with the primary-vortex core causing it to breakdown at 
an a n e lf °f attack as low as 18®. Such a critical angle 
of attack is substantially smaller than the critical angle 
of of vortex breakdown in the low-speed regime. 
Reference 21 contains extensive flow measurements for 
the 65° cropped delta wing with and without leading- 
edge extension (LEX). A complete reconstruction of the 
three-dimensional flow field at and behind the terminating 
shock was not possible experimentally. 

Computational simulations for transonic delta-wing 
flows have been developed on a very limited scale by 
using the Euler equations 20, 22 and the thin-layer Navier- 
Stokes equations 23 . The Euler-equations solutions were 
not capable of fully resolving the flow in the terminating 
shock region and the thin-layer Navier-Stokes-equations 
solutions did not address that region. In Ref. 24 by the 
present authors, the laminar, unsteady, compressible, full 
Navier-Stokes equations are integrated time accurately us- 
ing the implicit, upwind, flux-difference splitting, finite- 
volume scheme to study and construct the flow field 
structure of a transonic flow around a 65® sharp-edged, 
cropped-delta wing at 20® angle to attack, 0.85 Mach 
number and 3.23 xlO* Reynolds number. A fine O-H 
grid consisting of 125 x 85 x 84 points in the wrap-around, 
normal and axial directions, respectively, is used for the 
computational solution. A A-shock system, which con- 
sists of a ray shock under the primary vortex core and 
a transverse terminating shock, has been captured. Be- 
hind the terminating shock, the leading-edge vortex core 
breaks down into a two-bubble cell type. The terminat- 
ing shock and the vortex breakdown region behind it are 
time dependent and appear to be oscillatory. The flow 
field ahead of the terminating shock is steady and in- 
cludes a supersonic pocket which is surrounded by the 
ray shock and the terminating shock. The flow inside 
the pocket does not change due to changes in the flow 
downstream. This is consistent with the fact that the su- 
personic pocket along with the terminating shock do not 
allow disturbances to propagate upstream. These results 
have been validated using the available experimental data 
and they are in good agreement. This work gives a com- 
plete construction of the flow field over the wing surface 
and in particular the structure of the flow at the terminat- 
ing shock and behind it. 

In this paper, a parametric study is carried out to in- 
vestigate the effects of frecstream Mach number and an- 
gle of attack on the terminating shock and the leading- 
edge, primary-vortex breakdown for the same 65° sharp- 
edged, cropped delta wing. The computational investiga- 
tion uses the same equations, computational scheme and 
grid of Ref. 24. Keeping the Reynolds number fixed at 
3.23 x 10 6 , the Mach number is changed from 0.85 to 0.9 
while the angle of attack is fixed at 20°, and the angle 
of attack is changed from 20° to 24° while the Mach 
number is fixed at 0.85. 


HIGHLIGHTS OF FORMULATION 
AND COMPUTATIONAL SCHEME 

The conservative form of the dimensionless, unsteady, 
compressible, full Navier-Stokes equations is used for the 
formulation of the problem. The equations are written in 
terms of the time-independent, body-conformed coordi- 
nates and (Ref. 25). 

The implicit, upwind, flux-difference splitting, finite- 
volume scheme is used to solve the unsteady, compress- 
ible, full Navier-Stokes equations. The scheme uses the 
flux-difference splitting scheme of Roe which is based 
on the solution of the approximate one-dimensional, Rie- 
mann problem. In the Roe scheme, the inviscid flux dif- 
ference at the interface of computational cells is split into 
two parts; left and right flux differences. The splitting is 
accomplished according to the signs of the eigenvalues of 
the Roe averagcd-Jacobian matrix of the inviscid fluxes 
at the cell interface. The smooth flux limiter is used to 
e limina te oscillations at locations of large flow gradients. 
The viscous- and heat-flux terms are linearized in time 
and the cross-derivative terms are neglected in the im- 
plicit operator and retained in the explicit terms. The vis- 
cous terms are differenced using a second-order accurate 
central differencing. The resulting difference equation is 
approximately factored and is solved in three sweeps in 
the and directions. The computational scheme 
is coded in the computer program “FTNS3D” which is a 
modified version of the CFL3D-code. 

COMPUTATIONAL RESULTS 

A 65® swept-back, sharp-edged, cropped delta wing 
of zero thickness is considered for the computational so- 
lutions. The cropping ratio (tip length/root-chord length) 
is 0.15. An O-H grid of 125 x 85 x 84 in the wrap-around, 
normal and axial directions, respectively, is used. The 
computational domain extends two-chord length forward 
and five-chord length backward from the wing trailing 
edge. The radius of the computational domain is four- 
chord length. The minimum grid size normal to the wing 
surface is 5x10^ from the leading edge to the plane of 
symmetry. Figure 1 shows a three-dimensional shape of 
the grid and a cross-flow plane. 

Time-accurate integration of the laminar, unsteady, 
compressible, full Navier-Stokes equations has been car- 
ried out with At = 0.0002. Three flow conditions are 
used to study the effect of increasing the Mach number 
while the angle of attack is kept constant and the effect 
of increasing the angle of attack while the Mach num- 
ber is kept constant In all the three cases, the Reynolds 
number, R«, is 3.23x10* based on the root-chord length. 

Case I (Moo = 0.85, o = 20®) 

For this case, the frecstream Mach number. Moo. and 
angle of attack, o, are 0.85 and 20®, respectively. Figure 2 
shows a comparison of the computed, spanwise, surface- 
pressure coefficient (Cp) at different chord stations (x = 
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0.3, 0.6 and 0.8) with the experimental data of Erickson 
(R. = 3.23x10*) and Hartmann 17 (R« = 238x10* and 
4.57x10*). The computational results show the correct 
location and level of the suction-pressure peak corre- 
sponding to the primary vortex in comparison with the 
experimental data. They also show a smaller suction- 
pressure peak corresponding to the secondary vortex. The 
computational results are in fair to good agreement with 
the experimental data. For the chord station x = 0.9, the 
C p -curve shows a rapid increase in the pressure coeffi- 
cient (a decrease in the suction-pressure coefficient). For 
example, the suction-pressure-peak coefficient increases 
from a value of -1 .4 at x = 0.8 to a value of -1. 15 at x = 
0.9. Figure 3 shows the total-Mach contours and stream- 
lines at the chord stations of x = 0.60, 0.90 and 0.97. 
At x = 0.60, the Mach contours show an oblique shock 
beneath the primary vortex and a subsonic, separated re- 
gion to its right The streamlines show a secondary sepa- 
rated flow and the corresponding secondary vortex. This 
separation is due to the shock interaction with the sur- 
face boundary-layer flow and is also due to the adverse, 
spanwise pressure gradient created by the primary vor- 
tex. At x = 0.90, the shock beneath the primary vortex 
becomes weak and the primary-vortex size increases. At 
x = 0.97, the shock beneath the primary vortex disap- 
pears and the primary vortex diffuses and reduces to a 
repelling focus, as shown by the streamlines. The details 
of the flow structure at x = 0.90 and 0.97 in addition 
to the spanwise, pressure-distribution curve at x = 0.90 
clearly indicate that the primary vortex is experiencing a 
vortex breakdown due to a transverse shock (terminating 
shock) which is located between x = 0.80 and x = 0.90. 

Figure 4 shows the static pressure contours on the 
wing and symmetry planes. The contours clearly show 
the location, shape and strength of the terminating shock. 
A substantial supersonic pocket which is bounded by the 
terminating shock and the ray shocks (shocks beneath 
the primary-vortex cores) is observed on the wing plane. 
The terminating shock is located at x = 0.83 at the 
plane of symmetry, which is in good agreement with the 
experimental data 71 , where the shock is located at x = 0.84 
at the plane of symmetry. Figure 5 shows the position of 
ray lines from the wing vertex (which are marked by the 
letters A-H) and the static -pressure variation along these 
lines. The static -pressure curves give several points to 
generate the foot-print line of the terminating shock. The 
terminating shock is found to extend from the plane of 
symmetry to the wing leading edge. It reaches its highest 
strength at the location of the primary vortex (lines E-G). 
Figure 6 shows the total-Mach contours and streamlines 
on a vertical ray plane at the 0.68 spanwise location which 
passes through the vortex breakdown. Blow-ups of the 
velocity vectors and streamlines on this ray plane are also 
shown in Fig. 6. The streamlines conclusively show a 
two-bubble cell vortex breakdown. It is a typical three- 
dimensional vortex breakdown mode which consists of 
an attracting saddle point (front), a repelling saddle point 
(rear), an attracting focus (top), and a repelling focus 


(bottom). Such a breakdown mode is similar to the one 
which was captured for an isolated supersonic vortex in 
an unbounded domain in Refs. 26 and 27. The location 
of the attracting saddle point is at 0.97 along the ray line 
which corresponds to a location of 0.87 along the axial 
direction. The Mach contours show that the front surface 
of the vortex-breakdown bubbles is enclosed by a hemi- 
spherical shape-like shock surface. In Fig. 18, the details 
of the flow structure on the wing and symmetry planes 
are shown. 

Having established the flow structure of this case, the 
Mach number is increased to 0.9 while the angle of attack 
is kept fixed at 20°. 

Case II (Moo = 0.90, o = 20°) 

The results of this case are given in Figs. 7-1 1 and 
19. Figure 7 shows the computational spanwise, surface- 
pressure coefficient at different chord stations along with 
the experimental data of Erickson 21 . The computational 
results are in good agreement with the experimental data 
at x = 03 and 0.6. At x = 0.8, the computational re- 
sults underestimates the pressure coefficient of the exper- 
imental data. The locations of the primary and secondary 
vortex cores are in good agreement with those of the ex- 
perimental data. It is noticed that the levels of Cp for 
the present case are lower than those of Case I (Fig. 2). 
Again, the pressure level decreases rapidly at x = 0.90. 
Figure 8 shows the total-Mach contours and streamlines 
in cross-flow planes at x = 0.60, 0.90 and 0.97. The shock 
beneath the primary vortex is observed in the Figures at 
x = 0.60 and x = 0.90. For x = 0.90, the shock beneath 
the primary vortex is still strong in comparison with that 
of Case 1 (Fig. 3). At x = 0.97, the repelling focus is 
observed indicating that vortex breakdown has occurred. 
Figure 9 shows that the terminating shock in the cross- 
flow plane is located at x = 0.93 within the boundary- 
layer. which is in good comparison with the experimen- 
tally measured shock of Ref. 21, where it is located at x 
= 0.95. The static -pressure contours on the wing plane 
show that the terminating shock for Case II (Fig. 9) is 
closer to the trailing edge that of Case I (Fig. 4). It should 
be noted here that the terminating-shock location in the 
outer flow is ahead of its location in the boundary-layer 
flow. The static-pressure variations along the ray lines of 
Fig. 10 clearly show that the terminating-shock foot print 
is located between x = 0.925 and x = 0.95, and that it 
extends from the plane of symmetry to the wing leading 
edge. Figure 11 shows the Mach contours and stream- 
lines on a vertical ray plane passing through the vortex 
breakdown. It is noticed that the vortex breakdown shape 
is different from and smaller than that of Case I (Fig. 6). 
The attracting saddle point, attracting focus and repelling 
saddle point are clearly observed. The repelling focus 
is very small. This indicates that the terminating shock 
becomes smaller in strength than that of Case I. Figure 
19 shows the details of this flow case on the wing and 
symmetry planes. 
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It is that as the freestream Mach number 

increases slightly from 0.85 to 0.9, the terminating shock 
strength decreases and its location moves downstream 
from x = 0.84 to x ■ 0.93. Moreover, the surface pressure 
levels become smaller than those of Case I. 

Next, the Mach number is kept fixed at 0.85 and the 
angle of attack is increased to 24°. 

Case ID (Moo = 0.85, a = 24°) 

The results of this case are given in Figs. 12-17 and 
20. The computational surface-pressure result at x = 03 
(Fig. 12) is in good agreement with the experimental data 
of Erickson 21 . However, the computational results, at x 
= 0.6 and 0.8 are either overpredicting or underpredicting 
the experimental data. Figures 13, 14 and 15 show that 
the terminating shock moves upstream to x = 0.753 in the 
boundary-layer flow at the plane of symmetry. This is in 
good agreement with the experimental data of Ref. 21, 
where the shock is located at x = 0.75 in the boundary 
layer flow. The terminating-shock location in the outer 
flow is ahead of its location in the boundary layer. Figure 
16 shows that the vortex-breakdown region is larger than 
those of Cases I and II. Moreover, the attracting and 
repelling foci are smaller than those of Case I. Figure 20 
shows the details of this case on the wing and symmetry 
planes. 

Thus, it is seen that as the angle of attack increases 
from 20° to 24° while the Mach number is kept fixed 
at 0.85, the terminating shock moves upstream and the 
vortex-breakdown region becomes large. Moreover, the 
surface pressure levels become larger than those of Case I. 

The computational results show that the flow at the 
terminating shock and behind it is time dependent and it 
indicates oscillatory motion (The computations have not 
been carried out beyond t = 6.0 or 30,000 time steps 
with At = 0.0002). In Fig. 17, we show snapshots of the 
streamlines and their blow-ups on a ray plane passing 
through the vortex-breakdown region. The snapshots 
are shown at t = 4.22. 5.16 and 532. It is clearly 
seen that the vortex breakdown moves upstream showing 
different modes. In the same time, the terminating shock 
is also moving upstream and slows down to reverse its 
direction of motion. This is in complete agreement with 
the experimental observations of Bannik and Houtmann 16 . 

CONCLUDING REMARKS 

The laminar, unsteady, compressible, full Navier- 
Stokes equations are integrated time accurately using the 
implicit, upwind, flux-difTerence splitting finite-volume 
scheme to study the transonic flow field around a 65° 
sharp-edged, cropped delta wing. First, the flow field has 
been constructed for a Reynolds number of 3.23 xlO 6 , 
a Mach number of 0.85 and an angle of attack of 20® 
(Case I). A A-shock system consisting of a ray shock be- 
neath the primary vortex core and a transverse terminating 
shock has been captured. Behind the terminating shock. 


the leading-edge vortex core breaks down. Keeping the 
Reynolds number constant and the angle of attack fixed 
at 20®, the Mach number is increased to 0.90. The results 
of this (Case II) show that the terminating shock 
moves downstream and the vortex-breakdown region be- 
comes smaller than that of Case I. Keeping the Reynolds 
number constant and the Mach number fixed at 0.85, the 
angle of attack is increased to 20®. The results of this 
case (Case III) show that the terminating shock moves up- 
stream and the vortex-breakdown region becomes larger 
than that of Case I. The computational results arc in good 
agreement with the experimental data. However, it must 
be emphasized that the flow at the terminating shock and 
behind it is time dependent while the flow ahead of the 
terminating shock is steady. The present paper shows the 
structure of the flow field behind the terminating shock 
for the first time. 
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Fig. 1 Three-dimensional shape and cross-flow plane of a fine grid, 125x85x84. 
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Fig. 3 Total-Mach contours and streamlines in cross-flow planes; 
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Fig. 4 Static-pressure contours on the wing and symmetry planes, 

Moo = 0.85, a = 20°. 
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Fig. 6 Total-Mach contours, streamlines and velocity vectors on a ray plane 
passing through the vortex breakdown; Moo = 0.85, a — 20°. 
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Fig. 7 ‘ Comparison of the computed and experimental spanwise, surface-pressure 
coefficient at different chord stations; Moo - 0-90> a = 20°. 
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Fig. 12 Comparison of the computed and experimental spanwise, surface-pressure 
coefficient at different chord stations; Moo - 0.85, a = 24 . 
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Fig. 13 Total-Mach contours and streamlines in cross-flow planes; Moo - 0.85, a = 24°. 






PRESSURE CONTOURS ON THE WING SURFACE 


PRESS. CONTOURS ON THE PLANE OF SYMMTRY 


Fig. 14 Static-pressure contours on the wing and symmetry planes; 
Moo = 0.85, a = 24°. 
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Fig. 15 Ray lines on the wing surface and the static-pressure variation along them; 
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Fig. 18 Surface-pressure and Mach contours and particle trace on wing and 
symmetry planes; M 0 0 = 0.85, a = 20°. 






Fig. 19 Surface-pressure and Mach contours and particle trace on wing and 
symmetry planes; Moo ~ 0.90, = 20°. 
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Supersonic Vortex Breakdown on a Delta Wing 
M = 0.85, Re = 3,230,000 and AOA = 24 
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Fig. 20 Surface-pressure and Mach contours and particle trace on wing and 
symmetry planes; M oo = 0.85, a = 24°. 
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ABSTRACT 

The unsteady Euler equations and the Euler equa- 
tions of rigid-body dynamics, both written in the mov- 
ing frame of reference, are sequentially solved to simu- 
late the limit-cycle rock motion of slender delta wings. 
The governing equations of fluid flow and dynamics 
of the present multi-disciplinary problem are solved us- 
ing an implicit, approximately-factored, central-difference 
like, finite- volume scheme and a four-stage Runge-Kutta 
scheme, respectively. For the control of wing-rock 
motion, leading-edge flaps are forced to oscillate anti- 
symmetrically at prescribed frequency and amplitude 
which are tuned in order to suppress the rock motion. 
Since the computational grid deforms due to the leading- 
edge flaps motion, the grid is dynamically deformed using 
the Navier-displacement (ND) equations. Computational 
applications cover loc ally-conical and three-dimensional 
solutions for the wing-rock simulation and its control. 

INTRODUCTION 

The dynamic phenomenon of wing rock is character- 
ized by large-amplitude, high-frequency, rolling oscilla- 
tion with a limit-cycle amplitude. The rolling oscillation 
is self excited and it is triggered by vortex-flow asymme- 
try or vortex breakdown on highly swept delta wings at 
high angles of attack. The study of this phenomenon is 
vital for the dynamic stability and controllability of high 
performance aircraft during maneuvering and landing. 

The literature shows that several experimental 
investigations 1 ' 6 have been conducted to gain basic un- 
derstanding of the phenomenon. Nguyen, et al. 1 tested 
a flat-plate delta wing with 80° leading-edge sweep for 
forced-oscillation, rotary and free-to-roll tests. The free- 
to-roll tests showed that the wing exhibited a rock motion 
at angles of attack greater than 25°, and that the rock mo- 
tion reached the same limit-cycle response irrespective of 
the initial conditions. Levin and Katz 2 tested two delta 
wings with leading-edge sweeps of 76* and 80°. They 
found that only the wing with the 80° sweep would un- 
dergo a rock motion. Nelson and his co-workers 3 * 3 con- 
ducted a series of experimental studies to investigate the 
mechanisms responsible for wing rock on a delta wing 
with 80° leading-edge sweep. Their analysis revealed that 
the primary mechanism for the phenomenon was a time 
lag in the position of the vortices normal to the wing 
surface. Moreover, they concluded, through the analy- 
sis of separate contributions of the wing upper and lower 


surface-pressure distributions, that the upper surface pres- 
sure provides all of the instability and little damping in the 
roll moment and that the lower surface pressure provides 
the classical roll damping hysteresis. Morris and Ward 6 
conducted dynamic measurements in both a water tun- 
nel and a wind tunnel on a delta wing with leading-edge 
sweep of 80*. Their results showed that the measured 
hysteresis loops in the water tunnel were opposite in di- 
rection to those of the wind tunneL They concluded that 
the hysteresis direction does not play as decisive a role as 
previously thought in initiating and sustaining wing rock. 

Erickson 7,1 analyzed experimental data for aircraft 
configurations at high angles of attack in an attempt to 
reveal the flow processes which generate wing rock. He 
concluded that wing rock phenomenon for slender wings 
is caused by asymmetric-leading-edge vortices and that 
the vortex breakdown provides a limiter to the growth 
of wing-rock amplitude. He also identified another two 
mechanisms for limit-cycle oscillations in roll for ad- 
vanced aircraft. 

The literature review showed that numerical simula- 
tion of this phenomenon for tow speeds has recently been 
presented by Konstadinopoutos, et al.’. This has been 
followed by developments of analytical models to inves- 
tigate the parameters affecting this phenomenon. Nayfeh, 
et aL 1(Ml have presented two analytical models and Hsu 
and Lan 12 have presented one analytical modeL The im- 
proved analytical model of Nayfeh, et al. 11 proved to 
be superior in comparison with the Hsu and Lan model 
and more accurate than their first model of reference 10 . 
The model of reference 11 accurately fitted the rolling mo- 
ment coefficient, which was computed by a vortex-lattice 
method, using five terms which included the linear aero- 
dynamic damping and restoring moments and the nonlin- 
ear aerodynamic damping moments. With this model, it 
was shown on the phase plane that both the wing rock 
and wing-roll divergence were possible responses for the 
wing. Hsu and Lan’s model cannot predict wing-roll di- 
vergence. A serious question which can be raised regard- 
ing the work in references 9-12 is: how accurate the fluid 
dynamics solution is, using the vortex lattice method? 
Moreover, the fluid dynamics model limits its applica- 
bility to low-speed flows and to angles of attack below 
the critical value for vortex breakdown. Moreover, the 
vortex lattice model also cannot predict separated flows 
from smooth surfaces. 
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The first computational unsteady solution for the 
forced-rolling oscillation of a delta wing, which was 
based on the unsteady Euler equations, was presented 
by Kandil and Chuang 13 . The solution used the locally- 
conical flow assumption for supersonic flows in order to 
reduce the computational time by an order of magnitude 
as compared to that of the three-dimensional solutions. 
Forced-pitching oscillation of airfoils were also consid- 
ered in a later paper by Kandil and Chuang 14 . The first 
unsteady three-dimensional Euler solution for the forced- 
pitching oscillation of a delta wing was also presented 
by Kandil and Chuang 13 . The unsteady Navier-Stokes 
solutions were also used by Kandil and Chuang 16 for 
the forced-rolling oscillation of a delta wing under the 
locally-conical flow assumption. Batina 17 developed a 
conical Euler solver, which was based on the use of un- 
structured grids, and used it to solve for the flow around 
a delta wing undergoing forced-rolling oscillation under 
the locally-conical flow assumption. Later on, Lee and 
Batina' * extended the Euler solver to include a free-to- 
roll capability to solve for a freely rolling delta wing 
which exhibited wing rock. The solution was based on the 
locally-conical flow assumption. In Ref. 19, the present 
authors studied symmetric and anti-symmetric forced- 
rolling oscillations of the leading-edge flaps of a delta 
wing. A hinge is considered at the 75% location of the 
local half span and the leading-edge flaps are forced to 
oscillate both symmetrically and anti-symmetrically. The 
Navier-Stokes and Euler equations are used to solve the 
problem along with the Navier-displacement equation to 
account for the grid deformation due to the leading-edge 
flaps motion. In a later paper by the authors 20 , the effects 
of symmetric and anti-symmetric flaps oscillation with 
varying frequencies have been investigated for two flow 
conditions. With the aid of these studies, the authors 21,23 
studied the wing rock phenomenon as well as its ac- 
tive control using anti-symmetric tuned oscillations of the 
wing leading-edge flaps. The sequential solutions of un- 
steady Euler equations and the Navier-displacement equa- 
tions along with the Euler equation of rigid-body rolling 
motion were used to obtain the solutions for these prob- 
lems. The locally-conical flow assumption was also used 
throughout these solutions. Simulation of wing-rock and 
wing-divergence motions was presented by the authors 
for the three-dimensional flows in Ref. 23. 

In the present paper, the unsteady Euler equations and 
the Euler equations of rigid-body dynamics, both written 
in the moving frame of reference, are used to simulate 
the limit-cycle rock motion of slender delta wings. Con- 
trolling the wing-rock motion is achieved by using anti- 
symmetric forced-oscillation of the wing leading-edge 
flaps. For the active control of wing rock, the grid is 
dynamically deformed using the ND equations. 


FORMULATION 

The formulation of the problem consists of three sets 
of equations. The first set is the unsteady, compressible, 
Euler equations which are written relative to a moving 
frame of reference. This set is used to compute the 
flowfield for steady or unsteady flows. The second set is 
the unsteady, linearized, Navier-displacement equations 
which are used in the moving frame of reference to 
compute the grid displacements whenever the leading- 
edge flaps oscillate. If the leading-edge flaps do not 
oscillate, the ND equations are not used. The third set 
is the Euler equations of rigid-body motion for the wing 
only or for the wing and its flaps. This set is used to 
compute the wing motion for the wing-rock problem. It 
is solved in sequence with the first set For the control 
of wing-rock motion, this set is solved in sequence with 
the first and second sets. 

Unsteady Euler Equations 

Using the transformation equations from the space- 
fixed frame of reference to a moving frame of reference 
(Refs. 13-15), the non-dimensional, unsteady, compress- 
ible, Euler equations are transformed to the moving frame 
of reference. Such a transformation eliminates the mo- 
tion of the computational grid for rigid wings having 
time-dependent rigid-body motion. Since the flaps of the 
wings are allowed very small relative rigid-body motion 
per time step of the integration scheme, one must con- 
sider the computational grid as time-dependent whenever 
the grid is updated, and the grid speed in Eqs. (4) and 
(5) must be computed. Hence, the Euler equations are 
given by 


dQ 0F, _ * 

dt 


(i) 


where 


Q = flowfield vector 

q 1 

= J = -j\p . pu i » P«J> ptij, pc) (2) 

C 0) 


E m = inviscid flux 

= \lpUm,P*lUm +d l C‘P,P*lU m 

d£ m 

+ P«j Um + di£ m P, pUmh - -^-p]* (4) 

u m = <9^ tit + ^ (5) 
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S = source term due to rigid-body motion = — S 


+ (wxr) • a 0 + V 9 • (o t - wxV) + V ■ (u'xr) 


+ (wxr) • (<ixf)]}‘ 

(6) 

V = V. - V x = relative velocity 

(7) 

Vi = V„ + u>xr 

(8) 

o t = a„ + wxr + 2u>xV t + wx(wxf) 

(9) 


(10) 


over a short time range (t - 1,) where A, p and p are 
kept constants. This yields the equation 

-£vp* + Jgj- [|v(V • u) + V’u] 
da x 

= P -ft + <%(r) (13) 

In Eq. (12). we use R« m to refer to the mesh point 
Reynolds number which is different Grom the flow 
Reynolds number. This has been done in order to provide 
a limiter for the grid displacement to avoid grid distortion 
or overlapping, particularly in regions of high flow rever- 
sal Equation (13) is the vector form of the ND equations 
to be used for computing the grid-points displacement u 
subject to displacement boundary and initial conditions. 
The equation is a parabolic equation in time which is in- 
tegrated by using the alternating direction implicit (ADI) 
scheme. The constant C»(f) in Eq. (13) is computed from 
the preceding time-range integrations. 


. TP , V* V? 

p( 7 — 1) 2 2 
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The reference parameters for the dimensionless form 
of the equations are L } a x ,Lla x and p x for the length, 
velocity, time and density, respectively. Here, l is a 
reference length which is taken as the wing root-chord 
length. 

In Eqs. (l)-(ll), the indicial notation is used for con- 
venience. Hence the indices k , /, n and s are summation 
indices and m is a free index. The range of k, I, m, n, 
and s is 1-3 and dt 2 gf^. 

The term ^ represents the mth component of the 
grid velocity. It is set equal to zero when the grid is not 
being updated. In Eqs. (l)-(ll), p is the density, u n the 
relative fluid velocity component, V 0 and a, translation 
velocity and acceleration of the moving frame, V t and 
a, the transformation velocity and acceleration from the 
space-fixed to the moving frames of reference, w and u> 
the angular velocity and acceleration of the moving frame, 
r the fluid position vector, p the pressure, e and h the total 
energy and enthalpy per unit mass relative to the moving 
frame and 7 the gas index which is set equal to 1.4. 


Unsteady, Linearized Navier-Displacement 
Equations 

The details of the derivation of these equations are 
given by the authors in Ref. 20. The dimensionless form 
of these equations is given by 


— Vp + 


d 

R tm dt 


^V(V • u) + V 2 u 


= P 


0 J tj 

W 


( 12 ) 


where ti is the displacement vector of a grid point For 
each grid point (a fluid element), Eq. (12) is integrated 


Euler Equation of Rolling Rigid Wing With and 
Without Oscillating Leading-Edge Flaps: 

Figure 1 shows a sketch of a wing and its flaps which 
are undergoing rolling motions. The rolling motion of 
the flaps is anti-symmetric. The wing is fixed to an 
axile which rotates in bearings. The bearings damping 
coefficient is A. Torsional springs of stiffness k are 
assumed at the ends of the axle. The xyz axes which 
are fixed to the wing are assumed to coincide with the 
principal axes of inertia of the wing-flaps configuration. 
At section A-A, the wing half span is /j and the flap 
width is I 2 . The masses of the wing and each flap are mi 
and m2, respectively, and their respective mass-moment 
of inertias around their centers of mass are I cl and I ti . 
The generalized coordinates of the system are taken as 0i 
and $ 2 , which are measured from the horizontal position. 
If the aerodynamic moment of the wing and its flaps about 
the x-axis is C, and if one uses the Lagrangian dynamics 
for obtaining the governing equations of motion, one gets 
the following equation for the 8\ coordinate 

Cr — ^27„j — — mj/j/j coa 

+ TTljll fj02l Sin 021 

= ^7**1 + 27««2 — ~ m 2 hh cos 02i^ 0j 

— m 2 V20? sin 02i 

— 2012/1/2^1^11 sin 02i + A0i + kfi 1 (14) 

where 9j, = 0 2 - 0j , Ji„ and lu, are the mass moment 
of inertia of the wing and the flap, respectively, around 
the wing axis of rotation. If the angles 0i and 0u are 
assumed to be small then the linearized equation reduces 
to 

C r — (27,«2 — 2^/2 ~ 

= (7,., + 27.„-^-m2/i/,)0“, 

+ + k$\ 
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On the other hand, if the flaps are not deflected and the 
wing and its flaps roll as a rigid body, Eq. (15) becomes 

CV = i + A0 i + k$ i (16) 

where is the mass moment of inertial of the composite 
wing-flaps configuration without relative motion. 

Equation (16) governs the wing-rock problem while 
Eq. (15) governs the linearized control of wing-rock prob- 
lem by using a prescribed motion of the leading-edge 
flaps. 

COMPUTATIONAL SCHEMES 

The computational scheme used to solve Eqs. (1)- 
(11) is an implicit, approximately-factored, centrally- 
differenced, finite-volume scheme 13 ' 15 . Added second- 
order and fourth-order explicit dissipation terms are used 
in the difference equation on its right-hand side terms, 
which represent the explicit part of the scheme. The Ja- 
cobian matrices of the implicit operator on the left-hand 
side of the difference equation are centrally-differenced 
in space, and implicit second-order dissipation terms are 
added for the scheme stability. The left-hand side spa- 
tial operator is approximately factored and the difference 
equation is solved in three sweeps in the £*, £ 2 and £ 3 
directions, respectively. 

For the wing-rock problem, Eq. (16) is solved using 
a four-stage Runge-Kutta scheme. Starting from known 
initial conditions for 0 and 0, the equation is explicitly 
integrated in time in sequence with the fluid dynamics 
equations, Eqs. (1-11). Equation (16) is used to solve for 
0, 6 and 0 while Eqs. (1-11) are used to solve for C,. 
If the initial C T is nonzero, a case of asymmetric steady 
flow at initial conditions, the initial values of 0 and 0 are 
set equal to zero and the motion is initiated by the initial 
rolling moment 

For the control of the wing-rock problem using flaps 
oscillation, the motion of the flaps; 0 2J) hi and 0 2J are 
specified and Eq. (14) (nonlinear equation) or Eq. (15) 
(linearized equation) is used to solve for 0 it h and 0j. 
The fluid dynamics equations, Eqs. (l)-(ll), and the grid- 
deformation equation, Eq. (13), are sequentially used to 
solve for C,. 

COMPUTATIONAL APPLICATIONS 
AND DISCUSSION 

Simulation of Wing-Rock-Motion 
(Locally-Conical Flow) 

A delta wing of sweep-back angle of 80°, at an angle 
of attack of 35° and a Mach number of 1.4 is considered. 
The wing has an elliptic section with sharpened leading 
edges. The wing mass-moment of inertia about its x axis 
is 0.02, the bearing damping coefficient is 0.2 and the 
spring stiffness is 0.74. The unsteady Euler equations 


are solved for locafly-conical flows. The computational 
grid is of 64x64 x2 in the wrap around, normal and «i«i 
directions, respectively. For these flow conditions, the 
steady flow is asymmetric, and hence C, ? 0 at t * 0. 
Therefore, we set 0f = 0J = 0. The Euler equations of 
fluid flow and of rigid-body dynamics are sequentially 
integrated accurately in time with At » 0.0025. Figures 
2 and 3 show the results of this case. Figure 2 shows the 
time responses of 0j, C, and C» and the corresponding 
phase planes of h vz 0,, C r vz 0 t and C» vz 9 X . The time 
responses show the long time, t cr 7, it takes to build up 
the growing roll-angle response. The responses clearly 
show that the 9 X and Cr continuously increase in time 
with increasing frequencies. The limit-cycle response is 
reached at t ~ 21 which is clearly shown on the phase 
planes. The mean amplitude of 0i is -0.5°, its maximum 
is 40° and its minimum is -41*. Figure 3 shows snap 
shots of the surface-pressure coefficient and cross-flow 
velocity at the instants corresponding to points 1 and 2 
on Fig. 2. The strong asymmetric morion of the primary 
vortices are clearly seen. Also, the surface-pressure- 
coefficient response clearly shows the generation of the 
restoring rolling moment to the wing motion. 

Active Control of Wing Rock Using 
Leading-Edge Flaps Oscillation 

The next step is to control the wing rock response 
of the previous case. For this purpose a leading-edge 
flap hinge is assumed to be at the 76% location of 
the local-half-span length. The flaps motion is intro- 
duced at to * 13.02 when 0j ■ -4* and C, m 0.0. 
The flaps morion is anti-symmetric and is given by 
0ji(t) = BU sin kf(t - to) , where lq is the flap re- 
duced frequency. With the aid of the previous values of 
0i. C, and k of the wing (can be measured by sensors 
to feed back the leading-edge flaps motion), we chose 
02iau * -0.5* and kf ■ 6.7. Equation (15) for the wing- 
flaps motion is sequentially integrated accurately in time, 
with At ■ 0.0025, along with the Euler equations of fluid 
flow, and the ND equation is used for the grid deforma- 
tion. Figure 4 shows the time responses of 0j and Cr for 
the wing. It is clearly seen that 0 X response is damped 
within t - to * 13 with a mean value of 5*. However, 
the wing is still oscillating periodically around this mean 
position with a small amplitude. Next, the flaps motion 
is modified by dividing the amplitude 0 xlMAX by 1 + (t 
- to) so that it decays with time. Figure 5 shows the 
steady response of the wing at t - 30. The wing assumes 
an equilibrium position of 5° without any oscillation, lb 
check that this is a stable equilibrium position, the wing 
is disturbed at t « 40 with a small 0j. Figure 5 also shows 
the time responses of 0j and C, after the disturbance con- 
firming that the equilibrium position is stable. Figure 6 
shows the phase planes of the whole response history of 
0i and C,. Figures 7-9 show the same results as those of 
Figs. 4-6 when the same control is applied at U ■ 23.27, 
which is during the limit cycle response. 
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Simulation of Wing-Rock Motion (Three- 
Dimensional Flow) 

Next, we consider the three-dimensional-flow simula- 
tion of die wing-rock problem. 

A sharp-edged delta wing with a leading-edge sweep 
of 80° is considered for the computational applications. 
The angle of attack is set at 30° and the frees tream Mach 
number is chosen as 0.3 for low speed simulation. The 
wing mass-moment of inertia about its axis is 0.285, the 
bearings damping coefficient is 0.15 and the torsional 
springs stiffness is 0.74. The unsteady Euler equations 
are solved for the three-dimensional flows. The bound- 
ary of the computational domain consists of a hemispher- 
ical surface with it center at the wing trailing edge on 
its line of geometric symmetry. The hemispherical sur- 
face is connected to a cylindrical aftersurface with its 
axis coinciding with the wing axis. The hemispherical 
and cylindrical radii are two root-chord lengths and the 
downstream, circular exit boundary is at two root-chord 
lengths from the wing trailing edge. The grid consists of 
48x32x32 grid points in the wrap-around, normal and 
axial directions, respectively. The grid is generated in 
the crossflow planes using a modified Joukowski transfor- 
mation, which is applied at the grid-chord stations with 
exponential clustering at the wing surface. 

Since the steady flow solution is asymmetric, C T 
in Eq. (16) is of non-zero value and hence Eq. (16) is 
initially inhomogeneous. At t = 0, we set 6° = 6° = 0 
and release the wing with its initial M t value as the 
driving rolling moment. At t = At, Eq. (16) of the wing 
dynamics is integrated to obtain and hence and 
(At = 0.005). Then, Eqs. (1-11) of the fluid flow are 
integrated to obtain the components of the flowfield vector 
and hence p and C r . Next, t is increased to 2 At and the 
sequential integration of the dynamics equation and the 
fluid flow equations is repeated. The sequential solutions 
are repeated until the limit-cycle amplitude response is 
reached 

In Fig. 10, we show the roll angle, rolling-moment 
coefficient, C,, and normal-force coefficient, C*. versus 
time. Significant transient responses develop in the time 
range of t = 0 — * 22, wherein the amplitudes of the re- 
sponses increase and decrease. Thereafter, t > 22, the 
amplitudes of the responses continuously increase until 
t = 95. At t £ 95, the amplitudes and frequencies of 
the responses become periodic reaching the limit-cycle 
response. During the limit-cycle response, the maximum 
roll angle. 6 imix , is 10*. the minimum roll angle, 0inu»> 
is -11° and the period of oscillation is 3.53, which cor- 
responds to a frequency of 1.78. With At » 0.005, each 
cycle of oscillation in the limit-cycle response requires 
706 time steps. The shown responses, up to t = 140, 
required 28,000 time steps. 


Next, we consider one cycle of the limit-cycle 
response and analyze the roll angle, rolling-moment- 
coefficient and normal-force-coefficient responses to gain 
physical insight of the wing-rock phenomenon. For this 
purpose, we show in Fig. 11 0 lt C , and C. vz. f in 
the range of t = 135.19 — 138.72. This period of os- 
cillation is marked by the numbers 1, 2, 3, 4 and 5 in 
Fig. 11. In the first quarter of the cycle (1 — 2), the roll 
angle of the left side of the wing decreases from 0* — » 
-1 1° and the wing rolls in the clockwise (CW) direction, 
the rolling-moment coefficient increases and changes sign 
from -0.057 — ► 0.0 — ► + 0.023 and the normal-force co- 
efficient decreases and then increases from 2.68 — 2.65 

— 2.75. It is important to notice that the rolling moment 
changes its sign which means that the rolling moment 
during the first part of this quarter of the cycle is in the 
CW direction (the same direction as the motion) and in 
the second part of this quarter of the cycle is in the CCW 
direction (the opposite direction of the motion). Hence, 
the rolling moment increases the negative angle in the first 
part and then it limits the growth of the roll angle in the 
second part In the second quarter of the cycle (2 — 3) 
the roll angle increases from -11* — 0 and the wing rolls 
in the CCW direction, the rolling-moment coefficient in- 
creases and then decreases from 40.023 — 0.045 -* 0.04 
and the normal-force coefficients increases and then de- 
creases from 2.75 — ► 3.0 — 2.84. The rolling-moment 
coefficient is in the CCW direction (the same direction as 
the motion). In the third quarter of the cycle (3—4) the 
roll angle increases from 0 — 10° and the wing keeps its 
rolling motion in the CCW direction, the rolling-moment 
coefficient decreases and changes sign from 40.04 — 0 

— -0.038 and the normal-force coefficient decreases and 
then increases from 2.84 — 2.78 — 2.86. Again, it is no- 
ticed that the rolling moment changes its sign from CCW 
to CW directions and limits the roll angle growth. 

In Figs. 12 and 13, we show snapshots at points 2 and 
4, respectively; of the cross-flow-velocity vectors and the 
static-pressure contours at the chord stations of 0.54, 0.63 
and 0.79 and the surface-pressure coefficient at the chord 
stations of 0.54 and 0.63. In Fig. 12, the primary vortex 
on the right side is nearer to the upper wing surface than 
the one on the left side. Moreover, the primary vortex 
on the right is further away from the plane of geometric 
symmetry in comparison to the one on the left. The 
surface-pressure curves show large peaks on the right side 
and that the surface-pressure difference on the right side 
is larger than the one on the left side. This results into 
a CCW rolling moment at this maximum negative roll 
angle of -11*. In Fig. 13. the opposite process occurs; 
the surface-pressure difference on the left side is larger 
than die one on the right side and this results into a 
CW rolling moment at this maximum positive roll angle 
of 410*. These results are consistent with those of the 
experimental data of Refs. 3 and 4. 

In Fig. 14, we show the variations of the maximum 
static pressure of the vortex cores of the primary vortices 
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on the left and right sides versus the roll angle for the 
chord station of 0.54. The numbers on the figures cor- 
respond to those in Ftg. 11. Since the maximum static 
pressure of the core is proportional to the vortex-core 
strength, it is obviously seen that the primary vortex on 
the right side has a greater strength at point 2 as compared 
to that on the left side. The strength differential between 
the right and left vortices along with the locations of the 
vortex cores contributes substantially to the net total CCW 
rolling moment which limits the negative growth of the 
roll angle and reverses the wing motion. Similarly, it is 
concluded that the strength differential between the left 
and right vortices at point 4 substantially contributes to 
the net total CW rolling moment which limits the positive 
growth of the roll angle and reverses the wing motion. 

In Fig. 15, we split the rolling-moment coefficient 
into restoring and damping components similar to Kon- 
stadinopoulos, et al. 9 . First, the rolling-moment coeffi- 
cient C, is fitted using the following expansions in terms 
of 0 and 9: 

Cr = Oi9 + d 9 0 + Oj0 3 + Qt0 7 9 

+ dj0 7 9 + at0 3 + aj9 3 + a$6*0 

■f a 9 9 m 9 3 + Qiq0 7 9 3 -f* U| (12) 

The coefficients ai — a 12 are determined using a least- 
squares fit. A comparison of the original (-©-) and fitted 
(-*-) rolling-moment coefficients is shown in Fig. 15. 
Next, we split the fitted-rolling-moment coefficient into a 
restoring part, Af r , and a damping part. Mi, as follows: 

M r = + dj# 1 + 

+ (18) 

Mi = (a 2 + i i 4 0 2 + a t 9*)0 

+ (at + a 9 0 7 )0 3 + arf* (19) 

In Fig. 15, we also show M r and 9 versus time, and 
Mi and 0 versus time. Moreover, we show on these 
figures the numbers 1, 2, 3, 4 and 5 which correspond 
to the same numbers in Figs. 11 and 14. In the first 
quarter of the cycle (1—2), the roll angle 9 decreases 
from 0 — -11°, the restoring rolling moment becomes 
negative during the first part and positive during the 
second part and the damping rolling moment, which is 
negative at point 1, increases during the first part and 
becomes almost zero during the second part. It is very 
interesting to notice that M r and Mi are negative during 
the first part and hence they are in the same direction 
as the motion. During the second part, M r becomes 
positive reaching its maximum at point 2 when fl m ,. = 
-11° and hence it limits the angle growth. During the 
same second part. Mi becomes almost zero indicating a 
loss of damping rolling moment. In the second quarter 


of the cycle (2—3), M r stays almost constant during the 
first part and drops to zero in the second part when the 
roll angle becomes 0*. During the same second quarter. 
Mi continuously increases from 0 to a maximum positive 
value when the roll angle becomes 0. In the third quarter 
of the cycle (3-4), a similar interaction of 9, M, and 
Mi as that of die first quarter (1-2) occurs except with 
opposite signs. These conclusions are exactly similar 
to those of Ref. 9. Hence, the loss of damping rolling 
moment is responsible for the wing-rock motion. 

CONCLUDING REMARKS 

The multidisciplinary problem of wing-rock motion 
and its active control has been simulated using the un- 
steady, compressible, Euler equations; the Euler equa- 
tion of rigid-body dynamiocs and the ND equations for 
the grid deformation. The fluid flow Euler equations are 
solved using an implicit, approximately factored, central- 
difference, finite-volume scheme; rigid-body Euler equa- 
tion is solved using a four-stage, Runge-Kutta scheme and 
the ND equations are solved using an ADI scheme. Sim- 
ulation of the wing-rock problem is obtained for a delta 
wing which is mounted on an axle with torsional springs 
and the axle is free to rotate in bearings with viscous 
damping. The wing starts its motion under the effect of an 
initial rolling moment due to the initially asymmetric flow 
at zero roll angle and zero angular velocity. For the ac- 
tive control of wing-rock motion, a tuned anti-symmetric 
leading-edge flaps oscillation is used to achieve that pur- 
pose. Also, it has been shown that the hysteresis re- 
sponses of position and strength of the asymmetric right 
and left primary vortices are responsible for the wing rock 
motion. Moreover, it has also been shown that the loss 
of aerodynamic damping rolling moment at the zero an- 
gular velocity value is a main reason for the wing rock 
motion. These conclusions are consistent with the pre- 
vious findings of the experimental 3,4 and computational 9 
research work. 
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Fig. 10. Roll angle, roll-moment-coefficient and normal-force-coefficient response 
for wing-rock motion; delta wing, a = 30°, M w = 0.3, Ixx = 0.285, A = 
0.15, fc = 0.74, d* = 6* = 0. 
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Fig. 11. Time responses for wing-rock motion during the limit cycle response. 



Fig. 12. Snapshot at point 2 of crossflow velocity, static-pressure contours and 
surface pressure for wing-rock motion. 
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Fig. 13. Snapshot at point 4 of cross flow velocity, static-pressure contours and 
surface pressure for wing-rock motion. 
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Fig. 15. Splitting of rolling moment (C r ) into restoring rolling moment (M r ) and 

damping rolling moment (Mj) for wing-rock motion during the limit-cycle response. 











